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FOREWORD 


This  technical  report  is  an  adaptation  of  the  master's  thesis  of  the 
same  title  by  Lori  B.  Orenstein.  Ms.  Orenstein  was  enrolled  in  the 
Department  of  Mechanical  Engineering  of  The  University  of  Texas  at 
Austin  and  received  her  degree  in  August  1982. 

Since  the  printing  of  the  thesis  the  author  has  made  a  few  subtan- 
tive  revisions  that  are  included  in  this  report,  primarily  at  the  end  of 
Chapter  5  and  in  Chapter  6.  Correction  of  a  small  error  in  the  calcula¬ 
tion  of  dispersion  in  computer  program  LB01  (line  70;  see  Appendix  B) 
led  to  a  small  improvement  in  the  waveforms  of  the  computed  N  waves. 
Figures  3.2  and  5.5  (and  Table  5.2)  were  revised  accordingly.  Additional 
waveform  computations,  which  allowed  a  comparison  of  the  effects  of 
thermoviscous  absorption  by  itself,  relaxation  absorption  (and  dispersion) 
by  itself,  and  combined  thermoviscous  and  relaxation  absorption  (and 
dispersion),  led  to  two  new  figures.  Figs.  5.6  and  5.7,  and  a  discussion 
of  their  implications.  The  new  material  strengthens  the  main  conclusions 
in  Chapter  6. 

This  research  was  carried  out  at  Applied  Research  Laboratories  and 
was  supported  by  the  Office  of  Naval  Research  under  Contract 
N00014-75-C-0867.  Scientific  Officer  for  ONR  was  Dr.  Logan  E.  Hargrove. 

David  T.  Blackstock 
Supervisor 
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CHAPTER  1 

INTRODUCTION  AND  REVIEW  OF  THE  LITERATURE 

In  this  thesis  we  describe  an  experimental  study  of  the  rise 
time  of  spark  produced,  spherically  spreading  N  waves  of  finite  amplitude 
The  main  purpose  of  the  study  was  to  determine  the  effects  of  oxygen 
and/or  nitrogen  relaxation  on  the  rise  time. 

A.  Introduction 

The  idealized  pressure  signature  of  an  N  wave  is  shown  in 
Fig.  1(a).  The  wave  begins  with  a  condensation  headed  by  a  shock,  which 
raises  the  local  pressure  above  its  ambient  value.  The  amount  the  local 
pressure  exceeds  the  ambient  pressure  is  called  the  excess  pressure  and 
is  denoted  by  the  symbol  p.  Following  the  shock,  the  excess  pressure 
decays  linearly  to  a  rarefaction.  Finally,  a  tail  shock  returns  the 
local  pressure  to  its  ambient  value.  The  peak  value  of  the  excess  pres¬ 
sure  is  called  the  amplitude  and  is  denoted  by  the  symbol  P.  The  time 
elapsed  in  the  condensation  phase  is  called  the  half  duration  and  is 
denoted  by  the  symbol  T.  In  reality,  the  shock  portions  of  the  waveform 
have  a  finite  rise  time  tr>  A  more  realistic  pressure  signature  of  an 
N  wave  is  shown  in  Fig.  1(b). 

We  chose  to  study  only  the  condensation  because  for  our  N  waves 
that  can  be  generated  in  the  laboratory,  the  rarefaction  is  not  a  mirror 
image  of  the  condensation.  In  particular,  the  tail  shock  is  not  as  well 
formed  as  the  head  shock.  In  addition,  the  rarefaction  phase  is  more 

likely  to  be  contaminated  with  reflected  and/or  diffracted  signals. 
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FIGURE  1.1 

N  WAVE  PRESSURE  SIGNATURE 
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The  spark  produced  N  waves  resemble  the  sonic  boom  pressure 

signal  created  by  supersonic  aircraft.  Sonic  boom  waveforms  have  been 

observed  to  have  extremely  long  rise  times.  The  long  rise  times  have  been 

12  3  4 

attributed  to  inhomogeneity  of  the  atmosphere,  ’  turbulence,  ’  and 

5  6 

oxygen  and/or  nitrogen  relaxation.  *  In  our  spark  experiments,  which  are 

carried  out  indoors  in  a  laboratory,  atmospheric  inhomogeneity  and 

turbulence  are  of  course  not  present.  We  are  therefore  able  to  study  the 

effect  of  relaxation  by  itself.  There  is  another  important  difference 

between  our  spark  produced  N  waveforms  and  actual  sonic  boom  waveforms. 

Although  we  can  accurately  recreate  the  amplitude  of  a  sonic  boom,  the 

-4 

half  durations  we  observe  are  shorter  by  a  factor  10  . 

According  to  lossless,  linear  theory,  sound  waves  are 
propagated  with  constant  energy  and  uniform  wave  velocity.  However,  the 
acoustic  waveform  of  an  N  wave  is  in  general  altered  from  that  predicted 
by  lossless,  linear  theory  by  two  sets  of  physical  effects.  The  first 
set  consists  of  the  various  absorption  mechanisms.  These  include  the 
molecular  transport  processes  (viscosity,  heat  conduction,  and  species 
diffusion),  rotational  relaxation  of  nitrogen  and  oxygen  molecules 
(modeled  as  bulk  viscosity),  and  vibrational  relaxation  of  nitrogen  and 
oxygen  molecules  (caused  by  large  departures  from  equilibrium  in  some  of 
the  internal  degrees  of  freedom).  The  second  set  consists  of  the  non¬ 
linear  propagation  effects  described  in  finite  amplitude  theory.  These 
effects  cause  high  amplitude  sound  waves  to  distort.  There  is  a 
steepening  of  the  compression  phase  of  the  waveform  and  an  easing  of  the 
expansion  phase.  The  absorption  mechanisms  tend  to  counteract  the 
steepening  process  but  aid  the  easing  process.  If  the  nonlinear  effects 


4 

are  strong  enough,  shocks  form  but,  because  of  absorption,  they  have 
a  finite  rise  time. 

The  N  waves  used  in  this  study  were  created  by  sparks.  The 
sparks  were  made  by  applying  a  potential  across  a  capacitor  and  allowing 
it  to  discharge  through  an  air  gap  between  two  electrodes.  The  amplitude 
and  half  duration  of  the  N  wave  were  altered  by  changing  the  potential 
applied,  the  capacitance,  the  length  of  the  gap  between  the  electrodes, 
and  the  distance  between  the  spark  source  and  the  receiving  apparatus. 
Each  time  the  N  wave  was  altered,  we  measured  its  amplitude,  rise  time, 
and  half  duration.  The  results  were  compared  with  the  theoretical 
predictions. 

Most  theoretical  predictions  for  N  wave  rise  time  are  obtained 
from  solutions  for  a  steady  or  a  step  shock.  The  boundary  conditions 
for  a  step  shock  are  used,  rather  than  those  of  an  N  wave,  because 
the  equations  of  motion  are  much  more  difficult  to  solve  for  an  N  wave. 
However,  there  is  some  question  whether  it  is  valid  to  compare  our 
measured  N  wave  rise  times  with  rise  times  predicted  for  step  shocks. 

An  N  wave  only  approximates  a  step  shock  in  the  limit  as  the  half  dura¬ 
tion  goes  to  infinity.  If  the  half  duration  is  very  short,  as  in  our 
case,  a  step  shock  is  actually  a  poor  approximation  for  an  N  wave. 

Two  different  models  are  used  to  analytically  predict  the 
rise  time,  the  thermoviscous  model  and  the  relaxation  model.  In  the 
thermoviscous  model.  Burgers'  equation  is  chosen  to  describe  nonlinear 
propagation  in  a  quiet  thermoviscous  gas.  The  absorption  mechanisms 
considered  are  viscosity,  heat  conduction,  species  diffusion,  and  rota¬ 
tional  relaxation  of  nitrogen  and  oxygen  molecules.  In  some  cases 
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Stokes's  assumption  is  made  and  rotational  relaxation  is  ignored.  In  the 
relaxation  model,  a  "Burgers-like"  equation  is  chosen  to  describe  finite 
amplitude  propagation  in  a  relaxing  gas.  Because  the  form  of  this 
equation  is  known  only  for  a  mono-relaxing  gas,  the  most  dominant  mode  of 
relaxation  must  be  chosen  and,  since  the  relaxation  time  for  oxygen 
is  much  longer  than  that  for  nitrogen  or  carbon  dioxide  (except  when 
the  relative  humidity  is  less  than  10%),  oxygen  is  assumed  to  be  the 
dominant  mode  of  relaxation.  The  other  vibrational  modes  are  said  to  be 
in  equilibrium.  Although  thermoviscous  effects  may  be  included  in  the 
relaxation  model,  they  are  ignored  in  the  model  considered  here. 

Because  of  the  questionable  applicability  of  the  step  shock 
solutions  to  our  N  waves,  we  developed  a  computer  algorithm  which 
simulates  finite  amplitude  propagation  in  quiet,  but  real,  air.  The 
wave  propagation  algorithm  may  be  used  for  arbitrary  signals,  including 
our  N  waves.  The  ANSI  standard  SI. 26-1 978^  was  used  to  specify  the 
effects  of  viscosity,  heat  conduction,  rotational  relaxation,  vibrational 
relaxation  of  oxygen  and  nitrogen  molecules,  and  dispersion  due  to  oxygen 
relaxation.  The  theory  of  finite  amplitude  waves  was  used  to  determine 
the  nonlinear  propagation  effects.  Spherical  spreading  was  also  included 
in  the  algorithm.  Given  the  initial  waveform  and  the  propagation  dis¬ 
tance,  we  were  able  to  predict  the  local  waveform  and  also  the 
rise  time. 

B.  Literature  Review 

1 .  Theoretical  Literature 

The  theoretical  approach  to  nonlinear  propagation  in  the 
atmosphere  started  with  the  investigation  of  finite  amplitude  propagation 
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8  o 

in  a  lossless  fluid.  Stokes  (1848)  showed  that  Poisson's  solution 
implies  a  steepening  of  the  compression  parts  of  the  waveform  and  an 
easing  of  the  expansion  parts.  Using  an  arbitrary  curve  he  demonstrated 
that  the  peaks  travel  faster  than  the  zeroes  of  the  wave  and  the  troughs 
travel  slower.  As  time  increases  the  compression  portions  become  steeper. 
Eventually  a  discontinuity  is  formed,  at  which  time  the  equations  of  motion 
are  no  longer  valid.  Using  conservation  of  mass  and  momentum,  he  found 
conditions  under  which  the  discontinuity  can  be  propagated.  Rankine^ 
(1870)  found  that  unless  heat  conduction  from  particle  to  particle  is 
included,  there  is  a  loss  of  energy  in  the  propagation  of  a  finite  longi¬ 
tudinal  disturbance  that  is  not  accounted  for.  By  including  heat  conduc¬ 
tion  in  his  analysis  he  found  a  form  of  a  finite  longitudinal  wave  that 
propagates  without  changing  its  shape.  Rayleigh  ( 1 91 0 ) 1 1  found  that, 
although  heat  conduction  is  a  good  dissipation  mechanism  for  condensa¬ 
tion  shocks,  rarefaction  shocks  are  impossible.  The  reason  is  that 
heat  transfer  by  means  of  conduction  requires  that  heat  pass  from  a 
hotter  body  to  a  colder  body,  and  in  the  case  of  a  rarefaction  shock, 
heat  would  pass  from  a  colder  body  to  a  hotter  body.  Rayleigh  applied 

Rankine's  method  for  the  more  general  case,  where  viscosity  is  included. 

12 

Taylor  (1910)  used  different  arguments  to  come  to  the  same  conclusion 
as  Rayleigh,  namely,  that  only  compression  shocks  are  possible.  He 
obtained  an  approximate  solution  of  the  nonlinear  equations  of  motion 
for  a  viscous,  heat  conducting  gas.  His  solution  is  valid  for  weak 
shocks.  From  this  solution  arises  the  term  Taylor  shock  thickness.  The 
model  equation  for  waves  that  can  contain  shocks  in  a  thermoviscous  qas 
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where  t'=t-x/c0,  and 


c2u 

ox 


BUUt, 


Jt't' 


(1.1) 


(1.1a) 


The  small  signal  sound  speed  is  cQ,  B=(y+1)/2,  and  y  is  the  ratio 
of  specific  heats.  The  absorption  mechanisms  are  included  in  6.  Viscous 
absorption  is  represented  by  4v/3,  where  v=u/pq.  v  is  kinematic  vis¬ 
cosity,  u  is  the  shear  coefficient  of  viscosity,  and  pQ  is  the  ambient 
density.  Absorption  due  to  heat  conduction  is  represented  by  v(y-l)/Pr, 
where  Pr  is  the  Prandtl  number.  Rotational  relaxation  is  included  in  the 
term  v(uv/p),  where  py  is  the  bulk  viscosity. 

Burgers  (1948)  proposed  an  equation  of  the  general  form  of 
Eq.  (1.1).  Hopf  (1 950) 1 4  and  Cole  ( 1 951 ) 1 5  found  the  general  solution 
of  Burgers'  equation.  Lighthill  (T 956) ^ 6  showed  that  Burgers'  equation 
is  a  good  approximation  for  the  exact  wave  equation  for  a  thermoviscous 
medium. 

Polyakova,  Soluyan,  and  Khokhlov  (1962)^ 7  were  the  first  to 
derive  an  equation  for  finite  amplitude  propagation  in  a  relaxing  medium. 
The  equation  is  similar  in  form  to  Burgers'  equation. 


T 


+ 


-  mT  n 

- 

o 


t't' 


(1.2) 


where  t  is  the  relaxation  time  and  m  is  the  dispersion.  A  sound  wave 
passing  through  a  relaxing  fluid  may  cause  a  change  in  the  thermodynamic 
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equilibrium  of  the  fluid.  The  fluid  requires  a  finite  amount  of  time  (t) 
to  reestablish  equilibrium.  If  the  frequency  of  the  sound  wave  is  very 
high  the  fluid  cannot  maintain  perfect  thermodynamic  equilibrium.  In 
this  case  the  vibrational  energy  becomes  effectively  "frozen";  the 
fluid  behaves  as  if  there  were  no  vibrational  modes.  The  ratio  of  spe¬ 
cific  heat  y  reverts  to  its  classical  value,  y^=7/b.  The  sound  speed  has 
the  frozen  value  c^  because  of  its  direct  relation  to  Yf  If.  however, 
the  frequency  is  very  low,  the  changes  take  place  so  slowly  that  the 
vibrational  mode  is  always  in  equilibrium  and  the  ratio  of  specific  heats 
is  the  equilibrium  value 


7/2  +  c 

Ye  "  5/2  +  c  <  Yf 

where  c  is  the  vibrational  specific  heat  (cvl-b/R).  The  equilibrium 
sound  speed  (cQ)  is  directly  related  to  y6-  The  parameter  m  in  Eq.  (1.2) 
is  a  measure  of  the  dispersion  of  the  medium. 


Because  there  are  now  two  values  for  y,  there  should  be  two  values  for  6, 
i.e. , 


Ye  +  1  Yf  +  1 

Bg  ~  2  »  2 

Because  y6  and  Yf  are  almost  the  same,  and  6  is  the  coefficient  of  a  small 
term  (the  nonlinear  term),  Polyakova  et  al.  make  no  distinction  between 
the  two  8's.  Presumably,  they  use  6f 


4 
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Polyakova  et  al.1^  solved  Eq.  (1.2)  for  a  steady  step  shock. 
They  found  that  the  solution  is  continuous  for  only  a  small  range  of 
amplitudes.  The  restriction  on  amplitude  had  been  found  earlier  by 
Lighthill  (1956)J6  The  shock  velocity  is 


v 


c  +  4  u . 
o  2  jump 


» 


where  we  have  assumed  that  the  shock  is  traveling  into  a  quiet  fluid 
and  u.  is  the  magnitude  of  the  discontinuity  (or  j'ump)  in  particle 

J  Ulilp 

velocity.  Now,  as  long  as 

c0  <  v  <  cf  , 

the  shock  is  continuous.  A  shock  which  meets  that  restriction  is  called 
"fully  dispersed".  If  the  shock  velocity  is  greater  than  the  frozen  sound 
speed  it  is  called  "partly  dispersed".  Polyakova  et  al.  also  discussed 
the  problem  of  a  steady  step  shock  in  a  viscous,  relaxing  fluid.  No 
analytical  solution  for  the  profile  was  found.  Ockenden  and  Spence 

I  O 

(1969)  used  a  slightly  more  exact  set  of  equations  for  finite  amplitude 
propagation  in  a  relaxing  fluid.  They  made  a  distinction  between  6f  and 
6e-  The  equation  of  motion  simplifies  to  Eq.  (1.2)  when  Ye~Yf  Ockenden 
and  Spence  also  studied  the  steady  step  shock  in  a  viscous,  heat  conduct¬ 
ing,  and  relaxing  fluid,  but  found  no  solution  for  the  shock  profile. 
Hodgson  and  Johanneson  (1971  ,  1979)^’^  and  Lighthill  ( 1 956 ) 1  ®  assumed 
that  the  shock  waves  had  taken  shape  and  used  a  combination  of  the 
Rankine-Hugoniot  shock  relations,  the  equation  of  state,  and  the  relaxa¬ 
tion  equation.  They  presented  an  exact  solution  for  the  structure  of 
the  step  shock  with  no  further  assumptions  than  that  the  vibrational 
specific  heat  and  the  relaxation  frequency  are  constant.  The  actual 
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profile  was  found  using  a  numerical  method.  Morfey  (1979)  summarized 
Hodgson  and  Johanneson's  predictions  for  rise  time  and  maximum  allowable 
amplitude. 

All  of  the  aforementioned  work  on  relaxation  was  on  finite 

22 

amplitude  propagation  in  a  mono-relaxing  gas.  Clark  and  Rogers  (1964) 
took  up  the  problem  of  several  relaxing  modes.  They  found  that  if  the 
relaxation  time  of  one  mode  is  much  longer  than  that  of  the  others,  one 
may  assume  that  all  the  modes  that  have  short  relaxation  times  are  in 
equilibrium.  In  the  case  of  real  air,  the  relaxation  time  of  oxygen  is 
much  longer  than  that  of  the  other  modes  except  for  relative  humidities 
less  than  10%.  Since  relative  humidities  as  low  as  10%  are  rarely  found 
in  the  real  atmosphere,  we  assume  that  we  are  dealing  with  a  single 
relaxing  mode. 

All  of  the  solutions  mentioned  above  are  found  using  the 

boundary  conditions  for  a  step  shock.  Very  little  has  been  written  about 

the  solution  for  the  N  wave.  It  has  generally  been  assumed  that  the  step 

shock  is  a  good  approximation  for  the  sonic  boom  waveform.  Crighton  and 
23 

Scott  (1978)  solved  both  Burgers'  equation  and  the  equivalent  equation 
for  a  relaxing  gas  for  an  N  wave.  However,  the  solutions  are  not  in  a 
form  from  which  we  can  find  a  rise  time  prediction.  Honma,  Glass,  and 

O  A  pC 

Tsumita  (1981,  1982)  ’  have  solved  the  Navier-Stokes  equations 

numerically  including  vibrational  relaxation  of  oxygen  and  nitrogen 

for  explosions  of  a  pressurized  air  sphere.  They  compared  the  computed 

26 

waveforms  with  Hoi st- Jensen ' s  experimental  data  and  found  good 

agreement. 


2.  Experimental  Literature 


11 


There  have  been  very  few  laboratory  measurements  of  N  wave 

rise  time  in  air.  We  could  find  only  two  such  studies,  one  by  Wright27 

26 

and  the  other  by  Holst-Jensen. 

Wright  used  an  experimental  setup  similar  to  ours;  that  is,  he 
used  sparks  created  by  discharging  a  capacitor  across  two  electrodes. 

The  range  of  propagation  distances  used  was  10-200  cm.  Wright  reported 
data  runs  at  two  different  spark  energies.  The  N  wave  amplitude  at  a 
distance  of  20  cm  was  3.62  mbar  in  the  first  and  7.86  mbar  in  the  second 
run.  The  associated  half  durations  at  an  amplitude  of  0.65  mbar  and  an 
unspecified  distance  were  8.2  ysec  and  17.0  ysec.  The  range  of  rise  times 
he  obtained  was  0.5-2. 2  nsec.  Wright  compared  his  measured  rise  times 
with  those  predicted  using  the  thermoviscous  model  and  obtained  fair  agree¬ 
ment.  He  found,  however,  that  as  tho  amplitude  of  the  N  wave  decreased, 
the  measured  rise  time  fell  short  of  the  theoretical  prediction.  The 
departure  of  the  measured  rise  time  from  the  predicted  rise  time  began 
at  a  larger  amplitude  in  the  waveforms  with  shorter  half  durations. 

Holst-Jensen  used  both  sparks  and  exploding  wires  to  create 
N  waves.  He  was  able  to  obtain  N  waves  of  very  long  half  duration.  His 
experiment  was  conducted  in  a  very  large  room,  and  because  the  amplitude 
of  the  N  waves  decreases  very  rapidly  due  to  spherical  spreading,  he 
was  able  to  obtain  waveforms  with  small  amplitudes  and  long  half  dura¬ 
tions.  He  created  an  N  wave  with  a  large  amplitude  and  a  long  half 
duration  but  captured  his  first  waveform  for  analysis  at  4.08  m.  This 
process  allowed  the  amplitude  to  decrease  significantly  while  the  half 
durations  increased  slightly  (due  to  nonlinear  convection).  The  range 
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of  measurement  distance  was  4.08-21.60  m.  The  measured  amplitudes  were 
small  (<0.80  mbar).  The  N  waves  Holst-Jensen  recorded  were  more  like 
sonic  boom  waveforms  than  those  recorded  by  Wright  because  of  their  longer 
half  durations  and  smaller  amplitudes.  The  experimental  rise  times  he 
obtained  were  In  better  agreement  with  those  predicted  by  the  relaxation 
model  than  Wright's.  However,  like  Wright,  he  found  that  the  rise 
time  depends  on  the  half  duration  In  addition  to  the  amplitude. 

C.  Outline  of  the  Study 

Chapter  2  contains  a  description  of  the  analytical  models  used 
to  predict  rise  time.  The  computer  propagation  algorithm  is  discussed 
in  Chapter  3.  In  Chapter  4  we  describe  the  equipment  and  procedures  and 
In  Chapter  5  present  the  results.  Chapter  6  contains  our  conclusions  and 
recommendation  for  further  work.  Appendix  A  contains  derivation  of  the 
equation  of  state  for  a  relaxing  fluid.  Listings  of  the  driver  program 
and  subroutines  used  in  the  computer  propagation  algorithm  are  found 
in  Appendix  B. 


CHAPTER  2 


ANALYTICAL  PREDICTIONS 


A.  Introduction 

In  this  chapter  we  summarize  the  analytical  models  that  have 
been  used  to  predict  shock  rise  time.  Basically,  we  use  two  different 
models.  In  the  thermoviscous  model,  the  absorption  processes  that  affect 
the  shock  rise  time  are  determined  by  the  molecular  transport  processes 
(viscosity,  heat  conduction,  and  species  diffusion)  and  the  rotational 
relaxation  of  nitrogen  and  oxygen  molecules.  In  the  second  model, 
called  the  relaxation  model ,  the  absorption  is  determined  by  the  vibra¬ 
tional  relaxation  of  nitrogen  and  oxygen  molecules. 

The  general  theoretical  approach  found  in  the  literature,  for 
both  cases,  is  to  solve  the  conservation  equations  for  a  step  shock 
(a  shock  connecting  two  steady  flow  fields)  in  a  medium  having  the  appro¬ 
priate  dissipative  mechanism.  The  solution  thus  gives  the  profile  of 
a  steady  shock  wave.  Since  the  profile  is  known,  a  prediction  for  the 
rise  tine  can  be  found.  This  approach  can  be  expected  to  work  only  for 
shocks  that  are  almost  steady.  It  is  questionable  whether  the  shocks  in 
our  experimental  N  waves,  which  have  very  short  half  duration,  are  steady 
enough  to  approximate  step  shock.  An  alternative  approach  based  on  a 
computer  algorithm  that  may  be  used  for  waves  of  arbitrary  shape  is 
described  in  Chapter  3. 
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B .  Thermo v iscous  Mode 1 


As  noted  in  Chapter  1,  Burgers'  equation  is  a  reasonable 
approximation  of  the  exact  equations  of  motion  for  progressive  waves  in 
thermoviscous  fluids.  The  form  of  Burgers'  equation  that  is  useful  for 
experiments  in  which  the  time  waveform  of  the  sound  wave  is  measured  is 

coux  "  0uut'  =  2c~  ut't'  ’  (2,1) 

o 

where 


Viscosity,  heat  conduction,  species  diffusion,  and  rotational  relaxation 

are  included  in  the  "diffusion  coefficient"  6.  Several  authors  invoke 

Stokes's  assumption,  i.e.,  uv/u=0.  In  making  this  assumption  one  ignores 

the  rotational  relaxation  effects.  Stokes's  assumption  does  hold  for  the 

monatomic  gases,  but  not  for  air.  If  we  make  use  of  the  experimental 
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data  of  Greenspan  (1959),  we  find  that  for  dry  air  uv/u=0.6. 

Now  let  us  solve  Burgers'  equation  for  a  steady  shock.  There 
are  constant  conditions  (no  gradients)  far  ahead  of  and  far  behind  such 
a  shock.  Since  the  shock  is  steady,  it  travels  with  a  constant  velocity 
cQ.  The  shock  velocity  is  given  as16 


where  ufl  is  the  particle  velocity  ahead  of  the  shock  and  u^  is  the 
particle  velocity  behind  the  shock.  Choosing  v=cQ  requires  that  u^'-u^ 
If  all  the  parts  of  the  wave  travel  with  the  velocity  cQ,  then  u=u(t') 
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only.  Therefore  ux=0.  Our  boundary  conditions  are 


u  -*•  -u-j  as  t'  ■*  -® 

» 

(2.3a) 

u  ->■  +u^  as  t '  -v  +» 

9 

(2.3b) 

ut,  =0  at  t'  =  +« 

• 

(2.3c) 

Moreover,  we  choose  an  origin  such  that 

u  =  0  at  t'  =  0 

• 

(2.3d) 

Since  ux=0,  Eq.  (2.1)  becomes 

-euuf  =  ^  ut'f 

» 

which  may  be  integrated  once  to  obtain 


m 


ut 


u.,  =  C-, 


The  constant  may  be  determined  by  applying  the  boundary  conditions  of 
Eqs.  (2.3).  The  equation  becomes 


_i_  du_  =  2  _  2 
Bcq  dt1  U1  u 

This  equation  may  be  solved  by  separation  of  variables, 
substituted  the  result  is 


(2.4) 

If  Eq.  (2.3d)  is 


u  =  u.j  tanh( — t *)  .  (2.5) 

Equation  (2.5)  describes  the  profile  of  the  shock.  (A  good  approxima¬ 
tion  of  this  profile  is  given  by  the  sketch  labeled  D*100  in  Fig.  2.1.) 


/ 
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There  are  several  definitions  for  the  rise  time  in  popular  use 
(see  discussion  in  Chapter  4).  By  principle,  the  definition  we  use  is 


t 


r 


(2.6a) 


where  au  is  the  jump  in  the  magnitude  of  the  particle  velocity  and 

(du/dt’)  is  the  maximum  slope  of  the  head  shock.  In  this  case  au=2u,. 
'  max  1 

Note  that  a  closely  related  definition  for  rise  time  is 


t 


r 


AU 


/(**.) 

'midpoint 


(2.6b) 


where  (du/df  )m-j<jD0int  is  t^ie  s^ope  a  tan9ent  to  the  head  shock  at  its 
midpoint,  i.e.,  at  p=0.5P.  For  the  thermoviscous  model,  Eqs.  (2.6a)  and 
(2.6b)  are  the  same.  For  the  relaxation  model,  however,  the  two  equations 
are  different. 

To  find  the  rise  time  of  the  shock  in  a  thermoviscous  medium, 
we  first  note  that  the  hyperbolic  tangent  function  has  its  maximum  slope 
at  the  midpoint  (t'=0  in  Eq.  (2.5)),  where  u=0.  Equation  (2.4)  therefore 
yields 


m 


)  = 

/  du  \ 

'max 

\dt'  /, 

midpoint 


6cq(au) 

46 


(2.7a) 


Accordingly,  application  of  either  Eq.  (2.6a)  or  Eq.  (2.6b)  leads  to 


t 


r 


=  46 
6CQAu 
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If  we  want  to  find  t  in  terms  of  the  magnitude  of  the  jump  in  pressure, 
we  can  rewrite  au  using  the  impedance  relation  and  find  that 


_  V 

V  BAP 


(2.7b) 


In  general,  we  find  that  the  rise  time  is  inversely  proportional  to  Ap. 

In  addition,  the  rise  time  is  directly  proportional  to  the  diffusion 

coefficient  6  which  is  directly  proportional  to  the  kinematic  viscosity  v. 

Substituting  the  known  qualities  for  air  at  a  specific 

temperature,  we  can  find  an  expression  for  rise  time  in  terms  of  inverse 
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peak  pressure.  For  air  at  25°C, 


v  =  1.559  x  10‘5  m2/sec, 
pQ  =  1.18  kg/m3, 

8  =  1.2, 
y  =  1 .4,  and 
1/Pr  =  1.41. 

We  employ  two  predictions  based  on  the  thermoviscous  model.  In  the  first 
case,  TV1,  Stokes's  assumption  is  used  (uv/u=0).  In  the  second  case,  TV2, 
Greenspan's  data  is  used  (uv/u=0.6).  The  respective  rise  time  predictions 
are  as  follows: 


Case 

TV1 : 

5  *  2.958  x 

10"5 

tr  =  1.16/Ap 

Case 

TV2: 

6  =  3.893  x 

10"5 

tf  =  1 .53/Ap, 

where  Ap  is  given  in  mbar  and  tr  is  in  ysec.  In  other  words,  if 

Ap  *  1.0  mbar,  t  =  1.16  ysec  for  case  TV1  and  tp  =  1.53  usee  for 

case  TV2.  (A  plot  of  Ap*tr  for  Case  TV2  is  given  in  Fig.  2.2.)  The  rise 
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time  predictions  found  in  this  section  are  compared  with  our  experimental 


results  in  Chapter  5. 

C.  Relaxation  Model 

Our  analysis  for  a  relaxing  medium  basically  follows  that  of 
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Polyakova,  Soluyan,  and  Khokhlov  and  Rudenko  and  Soluyan.  They  com¬ 
bine  the  continuity  equation  and  the  momentum  equation  with  the  following 
equation  of  state  (see  Appendix  A  for  the  derivation). 


('  +  T  ar) 


0  +  t  jt) 

2  ,  ,  1  /aV\  .2 
V  +  2  r  2  )  p 

n0 

+  mrc 


2  dpi 
0  dt 


(2.8) 


where  m-^-Cg  is  the  dispersion.  After  making  the  appropriate 
approximations,  they  obtain  Eq.  (1.2),  which  is  repeated  here  for 
convenience. 


[ux  -  T  «v]t>  +  [UX  -  £  -f] 


Jl.  M 
2c  ut’t' 


A  steady  shock  solution  of  Eq.  (1.2)  is  easily  obtained.  The 
boundary  conditions  are  given  by  Eq.  (2.3a,  b,  c)  and  the  oriqin  is  set  by 
Eq.  (2.3d).  Again,  by  choosing  the  shock  velocity  to  be  cQ,  we  restrict 
u  to  be  a  function  of  t'  only.  As  a  result  the  ux  terms  drop  out  and  we 
obtain 


=  0 


Rearranging  and  integrating  once,  we  find  that 


uut, 


me 

IT  V 


(2.9) 


« 


t 


4 
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Substituting  the  boundary  conditions  into  Eq.  (2.9),  we  find  that 
2 

C^=u^/2t.  Equation  (2.9)  can  now  be  written  as 


Expanding  the  left-hand  side  and  integrating  once  we  obtain 


(2.10) 


-,n(uru2) ♦  ^  =  t  +  4  • 

If  we  invoke  Eq.  (2.3d),  C2*-£n(u^).  Substituting  this  result  gives  us 


or 


where 


t  *  + ^  « 


me 


(^) 


f-  = 


U-, 


\u1  -  u 


/U1  +  u\ 

+  Dan(  — - - 

VI  '  7 


(2.11) 


D  = 


me 
_ o 

6AU 


(2.11a) 


Equation  (2.11)  may  be  written  as 
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t'  =  tn 

l 

♦  »/u,r  i 

—  Jen 

T 

.1 

Equation  (2.12)  describes  the  profile  of  the  shock.  The 
profile  has  a  different  shape  depending  on  the  value  of  the  parameter  D, 
which  is  a  ratio  of  the  parameters  representing  dispersion  ar»d  non¬ 
linearity.  In  Fig.  2.1  the  shock  profile  is  shown  for  various  values  of 
D.  In  computing  these  profiles,  we  did  not  locate  the  zero  of  the  wave 
at  t'=0.  Another  starting  point  was  chosen.  As  a  result,  the  t'  origin 
varies  from  plot  to  plot.  If  D>>1 ,  the  nonlinear  effects  are  weak  and 
the  profile  becomes  a  hyperbolic  tangent  function  (i.e.,  the  shock  pro¬ 
file  is  perfectly  skew  symmetric).  The  profile  is  the  same  as  that  of 
a  shock  in  a  thermoviscous  medium  (see  Eq.  (2.5)).  As  D  decreases 
but  is  still  greater  than  1,  the  shape  of  the  profile  begins  to  lose  its 
perfect  skew  symmetric  shape.  The  shock  profiles  for  D»1  and  D>1  are 
called  fully  dispersed.  At  D=1 ,  there  is  a  cusp  at  u=-u-j  .  For  D<1  the 
profile  becomes  multivalued.  Since  a  multivalued  waveform  is  physically 
impossible,  we  expect  that  the  profile  becomes  discontinuous.  In  this 
case  the  shock  becomes  frozen  and  is  called  partly  dispersed.  We  can 
use  Eq.  (2.12)  to  predict  the  rise  time  for  fully  dispersed  shocks  only. 

If  thermoviscous  effects  are  included,  the  discontinuous  section  in  a 
partly  dispersed  shock  becomes  continuous.  However,  we  have  not  attempted 
a  solution  for  both  effects. 

The  expression  for  D  is  found  in  Eq.  (2.11a).  If  we  invoke 
the  impedance  relation,  we  can  find  D  in  terms  of  the  pressure  jump  Ap, 
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D  = 


mPQC0 

6Ap 


We  can  determine  the  maximum  allowable  peak  pressure  for  which  we  can 
predict  the  rise  time  by  setting  D=1  and  solving  for  Ap, 


AP„ 


m  r2 

mp  c 
o  o 


rmax  a 

Since  the  dispersion,  density,  and  sound  speed  change  with  temperature, 

we  need  to  choose  a  particular  temperature  for  the  calculation.  In  our 

experiment  the  temperature  was  fairly  constant  at  25°C.  The  value  for  m 

5  6 

is  calculated  by  using  the  harmonic  oscillator  equation  ’  to  calculate 
the  vibrational  specific  heat  of  oxygen.  From  this  result  we  fin  y e- 
The  dispersion  is  defined  as 


cr  -  c 
_  _  f  o 
m  -  — 2 — 


or 


'if  '  yc 


The  value  for  y^  is  always  1.4.  For  25°C 


-4 


ye  =  1.3989 
m  =  7.39  x  10 
pq  =  1.18  kg/rn^ 


cQ  =  345.94  m/sec,  and 


6  =  1.205. 

At  this  temperature,  the  maximum  pressure  ,iump  is  (oxygen  relaxation  only) 

Ap  =  0.866  mbar. 

Kmax 
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Above  this  limit  the  shocks  are  partly  dispersed.  We  cannot  compare  any 
data  for  pressure  jumps  above  0.866  mbar  with  the  predictions  from  this 
model . 

Recall  the  discussion  of  the  definition  of  rise  time  on  p.  16. 

If  the  rise  time  based  on  the  slope  at  the  midpoint  is  calculated,  we 
use 


midpoint 


2t(u+DUi ) 


t'=0 


AU 

4x0 


and  use  of  Eq.  (2.6b)  leads  to 


t 


r 


4tD 


4rmc 
_ o 

6AU 


(2.13a) 


(2.13b) 


If  we  want  to  find  t  in  terms  of  Ap,  the  magnitude  of  the  jump  in 
pressure,  we  can  rewrite  Au  using  the  impedance  relation  and  find  that 


t 


r 


4Tmpoco 

eAp 


(2.14) 


On  the  other  hand,  the  rise  time  calculated  by  means  of  Eq.  (2.6a) 
requires  the  value  of  the  maximum  slope,  which  is  found  to  be 


(#L  ■  *  (D  •  ■  <2-,5a) 

Substitution  of  this  value  in  Eq.  (2.6a)  leads  to  (after  some  manipula¬ 
tion) 


(2.15b) 


Now  consider  the  limiting  cases.  The  limit  as  D-*»  is 
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t  =  4tD 
r 

in  agreement  with  Eq.  (2.13b),  as  one  would  expect  (see  the  plot  for 
D=100  in  Fig.  2.1).  At  the  other  extreme,  D=1 ,  the  rise  time  is 

tr  =  2t 

Thus,  when  D=1 ,  the  value  of  the  rise  time  based  on  the  maximum  slope  is 
half  that  based  on  the  slope  at  the  midpoint  of  the  shock. 

The  choice  between  Eqs.  (2.13b)  and  (2.15b)  was  eventually  made 
on  the  basis  of  experimental  data.  For  the  N  waves  we  measured,  we  found 
that  the  maximum  slope  of  the  head  shock  occurred  at  the  shock  midpoint. 
It  therefore  seemed  reasonable  to  use  Eq.  (2.13b)  as  the  expression  for 
the  relaxation  model  and  Eq.  (2.6b)  as  the  definition  on  which  the 
experimental  measurement  is  based. 

In  Fig.  2.2  we  plot  the  rise  time  times  Ap  versus  the  relative 
humidity  for  several  temperatures  (based  on  Eq.  (2.14)).  For  comparison 
the  curve  TV2  is  included.  The  rise  time  is  heavily  dependent  on  the 
humidity  and  temperature.  The  dependence  is  mainly  a  result  of  the 
dependence  of  the  relaxation  time  on  humidity  and  temperature.  For 
example,  at  35%  relative  humidity,  15°C,  the  relaxation  time 
t  =  9.38  usee.  At  100%  humidity,  30°C,  the  relaxation  time 
t  =  0.931  ysec.  The  dispersion  m,  the  density,  and  the  equilibrium 
sound  speed  have  a  slight  dependence  on  temperature.  In  order  to  compare 
the  predictions  from  the  model  with  the  experimental  data  we  choose  the 
relative  humidity  and  temperature  and  find  Ap*tr>  For  example,  at 
60%  relative  humidity,  25°C, 

Ap  •  tr  =  9.638  , 
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where  Ap  is  in  mbar  and  t  is  in  n- sec.  Under  these  conditions  the  rise 
time  is=six  times  longer  in  a  relaxing  medium  than  in  a  thermoviscous 
medium.  The  appropriate  representation  for  t  is  used  for  comparison  with 
our  experimental  data  in  Chapter  5. 

D.  Discussion 

Before  comparing  our  experimental  rise  times  with  those 
calculated  using  the  models  described  in  this  chapter,  we  must  decide 
whether  the  models  fit  our  experiment.  First,  the  rise  time  predictions 
are  for  plane  waves,  whereas  our  experimental  data  are  for  spherical 
waves.  We  assume  that  the  plane  wave  formulas  are  valid  if  the  amplitude 
of  the  N  wave  is  changing  slowly  enough  that  the  shock  can  adjust  rapidly 
enough  to  follow  the  change.  This  assumption  may  be  poor  close  to  the 
source,  but  it  is  valid  at  great  distances.  Second,  and  more  important, 
the  theoretical  solutions  are  for  a  step  shock.  This  assumption  is  valid 
for  sonic  booms  because  the  sonic  boom  waveform  has  an  extremely  long 
half  duration  (~0.2  sec).  The  slope  behind  the  head  shock  is  very,  very 
gentle.  Our  N  waves,  however,  have  very  short  half  durations 
(-0.02  msec).  The  slope  behind  the  head  shock  is  very  steep,  and  may 
therefore  affect  the  rise  time. 

To  obtain  better  solutions  we  would  need  to  derive  the  rise 

time  predictions  directly  for  an  N  wave.  Since  the  shock  profile  is  not 

steady  (the  amplitude  decreases  linearly  behind  the  head  shock),  the 

simplification  ux=0  cannot  be  used.  Instead,  the  full  equation  (Eqs.  (2.1) 

or  (1.2))  must  be  solved.  Both  Burgers'  equation  and  Eq.  (1.2)  have  been 
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solved  for  an  N  wave.  Howeve. ,  in  neither  case  was  a  rise  time  predic¬ 
tion  made.  As  an  alternative  to  attempting  an  analytical  solution,  we 


used  a  computer  propagation  algorithm  similar  to  those  described  by 
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Pestorius,  Anderson,  Theobald,  and  Essert.  The  algorithm  is 

described  in  the  next  chapter. 


CHAPTER  3 


WAVE  PROPAGATION  ALGORITHM 


A.  Introduction 

In  this  chapter  we  present  our  alternative  to  the  analytical 
models,  the  wave  propagation  algorithm.  The  algorithm  is  based  on  the 
assumption  that  over  a  small  enough  propagation  distance  the  effects  of 
nonlinear  distortion  and  absorption  are  independent  and  may  therefore  be 
computed  separately  and  then  added.  The  algorithm  is  applicable  to 
arbitrary  waveforms,  even  ones  containing  discontinuities;  therefore 
it  is  ideal  for  dealing  with  our  N  waves.  It  is  designed  to  calculate 
the  shape  of  the  waveform  as  a  function  of  distance.  The  computer 
language  used  in  the  algorithm  is  FORTRAN.  The  programs  were  run  on  a 
COC  Cyber  Series  170  computer. 

The  theoretical  basis  of  the  algorithm  is  presented  in  Section  B 
below.  The  calculations  that  constitute  the  "absorption  step"  and  the 
"nonlinear  distortion  step"  are  given.  In  Section  C  we  examine  the 
algorithm  in  general.  The  driver  program  is  discussed,  part  by  part,  and 
the  subroutines  explained.  Two  sample  computed  waveforms  are  shown  in 
Section  D. 

B.  Theory 

1 .  Absorption 

Consider  a  single  frequency  component  of  an  arbitrary  spherical 

wave.  If  the  (complex)  amplitude  of  the  component  is  UQ  at  a  distance 

r  ,  the  amplitude  U  at  a  distance  r  is,  in  the  absence  of  nonlinear 
o 
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distortion,  given  by 


29 


-  roUn 

U  =  -p-  e  0 


The  factor  e 


-jk(r-r  ) 


does  not  appear  because  the  wave  is  represented  in 


the  time  domain  in  terms  of  the  retarded  time  t'=t-(r-r  )/c  .  If 

o  o 

r=rQ+Ar,  where  Ar  is  the  propagation  distance  step. 


U  = 


ro  +  Ar  o 


U  e 


-c(r-r  ) 


where 

C  =  a  -  jn  » 

a  is  the  attenuation,  and  n  is  the  dispersion  parameter. 


■  -e  -  i) 


(3.1) 


(3.1a) 


(3.1b) 


2/  3 

For  a  thermovi scous  medium,  u=6u>/2 cq  and  n=0.  For  a  relaxing  medium  the 
absorption  is 

2 


mxu) 

°  T~Z 


2c  (1+0)  T  ) 
o 


where  m=(c2-c2)/c2  and  the  dispersion  associated  with  it  is 


n  =  oixa 

For  the  real  atmosphere,  a  is  given  by  the  ANSI  standard,  which  is  dis¬ 
cussed  in  the  section  on  application  of  absorption.  The  dispersion  is 

n=o>(To)n  .  The  application  of  Eq.  (3.1)  is  called  the  absorption  step. 
°2 

It  includes  spherical  spreading,  as  well  as  absorption  and  dispersion. 


2.  Nonlinear  Distortion 


The  distortion  of  the  waveform  is  caused  by  a  relative  time 
shift  of  the  wavelets.  A  wavelet  is  a  point  on  the  waveform.  The  time 
shift  is  due  to  the  dependence  of  the  wavelet  propagation  speed  dx/dt  on 
the  wavelet  amplitude  u,  in  particular 

,  (3.2) 

where  6=( y+1 )/2  and  y  is  the  ratio  of  specific  heats.  Given  a  time 


ft  '  co  *  Bu 


waveform,  at  a  distance  r=rQ,  consider  the  typical  wavelet  u^ ,  whose  time 
coordinate  is  t^.  The  value  of  t.  at  r=rQ+Ar  is 

t  =  t  +  _AT  . 

S'  S  dx/dt 
new  old 


=  t. 


Ar 


*’old 


co  +  6ui 


(3.3) 


or,  if  the  binomial  expansion  is  used. 


Ar  /  eui 

t.  =  t.  +  ~  ( 1  -  ~ ■  +  higher  order  terms 
’new  ’old  o  \  o 


) 


.  „  Su.Ar 

.  Ar  i 

~  S  c - 2 

’old  o  c„ 


(3.4) 


When  we  compute  the  new  time  waveform,  we  shift  the  entire  waveform  by 
an  amount  Ar/cQ.  That  is,  we  express  u  in  terms  of  the  retarded  time 
t'=t-(Ar/co) .  Equation  (3.4)  now  becomes 


f  =  t  -  sy 
S  S  2 

new  old  cQ 


(3.5) 


Thus,  after  a  distance  step  Ar,  the  wavelet  appears  at  time  tj  instead 

new 

of  t.  .  We  call  the  application  of  Eq.  (3.5)  the  distortion  step.  We 
’old 

use  a  plane  wave  relation,  Eq.  (3.2),  to  calculate  the  distortion  for  a 
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spherical  wave.  The  justification  for  this  is  that  although  the 
amplitude  u.  is  fixed  during  the  nth  distortion  step,  it  drops  because  of 
absorption  and  spherical  spreading  during  the  N+lth  absorption  step. 

During  the  n+lth  distortion  step,  therefore,  u..  is  smaller  by  the  appro¬ 
priate  amount. 

C.  Wave  Propagation  Algorithm 

In  general,  the  wave  propagation  algorithm  involves  the  imple¬ 
mentation  of  Eqs.  (3.1)  and  (3.5)  alternately  in  small  steps  (&r).  The 

initial  waveform  is  input  at  r=r  .  The  waveform  is  transformed  from  the 

o 

time  domain  to  the  frequency  domain,  via  a  fast  Fourier  transform  (FFT), 
and  absorption  is  applied  (Eq.  (3.1))  over  the  first  step.  The  waveform 
is  then  transformed  back  to  the  time  domain  by  an  inverse  FFt|(FFT)_1J 
and  nonlinear  distortion  is  applied  (Eq.  (3.5))  over  the  first  step.  The 
propagation  distance  r  becomes  (rQ+Ar),  and  the  process  is  repeated  for 
the  next  step.  The  loop  continues  until  r  equals  a  specified  propagation 
distance. 

Several  important  details  have  been  left  out  of  the  discussion 

above.  First,  it  is  important  that  we  take  small  steps  to  ensure  that  no 

shocks  form  during  a  nonlinear  distortion  step.  We  apply  the  absorption 

step  before  the  distortion  step  to  take  care  of  the  possibility  that  the 

initial  waveform  contains  a  true  discontinuity.  Open  air  absorption  is 

sufficient  to  prevent  a  multivalued  waveform,  assuming  the  step  is  small. 

We  ensure  that  the  step  size  is  small  by  calculating  the  shock  formation 

distance  after  each  absorption-distortion  step.  We  then  use  this  calcu- 
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lation  to  select  the  next  value  of  Ar.  Pestorius  found  that  a  step 
size  of  one- tenth  the  shock  formation  distance  is  sufficiently  small 


32 


to  prevent  shocks  from  forming  in  the  subsequent  steps.  Although 
Pestorius  looked  at  plane  waves,  not  spherical  waves,  we  assume  that  his 
findings  are  applicable  to  our  problem.  We  also  improve  the  efficiency 
of  our  program  in  recalculating  the  step  size,  by  taking  advantage  of  the 
fact  that  Ar  can  be  increased  as  the  wave  gets  weaker.  As  a  result,  we 
can  reduce  the  number  of  trips  back  and  forth  from  the  time  domain  to  the 
frequency  domain.  The  second  omission  is  that  after  thr  distortion  step 
the  time  increments  between  the  points  in  the  waveform  are  no  longer 
equally  spaced.  The  FFT  routine  we  use  requires  that  the  time  increments 
be  equally  spaced.  We  find  equally  spaced  time  increments  using  the 
subroutine  RESAMP. 

The  wave  propagation  algorithm  is  implemented  by  a  driver 
program,  LB01 .  A  listing  of  the  program  is  found  in  Appendix  B.  A 
flowchart  of  the  algorithm  is  shown  in  Fig.  3.1. 

1 .  Driver  Program 

a.  Input  Scheme 

The  input  to  the  program  is  the  time  waveform  of  the  N  wave 
received  at  the  initial  position  r  .  A  data  file  is  created  by  reading 
the  voltage  levels  from  the  digital  oscilloscope  in  0.5  ysec  time  incre¬ 
ments.  A  buffer  of  zeroes  is  included  on  either  side  of  the  data  to 

29 

prevent  an  end  point  problem  similar  to  that  encountered  by  Pestorius. 

The  number  of  data  points  (N)  in  the  waveform  is  the  total  duration 
(including  zeroes)  divided  by  the  time  increment  (At). 

At  the  end  of  the  data  file  we  list  several  constants  used  in 
the  program.  Other  constants  are  input  on  the  terminal.  Some  of  the 
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important  parameters  and  constants  are  as  follows. 
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P  the  real  one-dimensional  amplitude  array  containing  N 
amplitude  coordinates  (u^  of  the  waveform  (in  mV). 

A  the  complex  amplitude  array  containing  the  frequency 
components  of  the  waveform  (used  only  in  the  frequency 
domain) . 

T  the  real  one-dimensional  time  array  containing  N  time 
coordinates  (t. )  of  the  waveform  (in  seconds). 

X  the  real  one-dimensional  array  describing  the  source- 
receiver  distances  (including  the  input  distance)  at 
which  we  want  an  output. 

NP  the  number  of  distances  in  array  X. 

N  the  number  of  points  in  the  waveform. 

DELT  the  time  increment  between  points. 

I OPT  the  plotting  option. 

FCON  the  conversion  factor  for  converting  the  P  array  to  a 
particle  velocity  array. 

PNORM  a  normalizing  factor  used  in  the  plot  routine. 

RN  the  propagation  distance. 

SSI  the  incremental  step  size  (&r)  used  to  cover  the 

distances  X(I). 

b.  Step  Size  Calculation 

The  step  size  is  calculated  for  each  pass  of  the  absorption  and 
nonlinear  distortion.  This  allows  us  to  reduce  the  calculation  time  by 
taking  advantage  of  the  slowing  down  of  the  distortion  caused  by  spherical 
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spreading.  The  step  size  is  calculated  using  a  relation  for  the 
distortion  variable. 


_  z 

a  -  —  » 

z 

where  z=r0$,n(r/r0)  and  z  is  the  value  of  z  when  shock  formation  occurs 
(a=l).  The  shock  forms  at  the  place  in  the  waveform  where  the  waveform 
has  its  maximum  slope.  We  find  that 


z  = 


elf7) 


max 


where  f‘  is  the  slope.  Combining  the  information  above  we  find  that 

r 


max 


infc)  ' 


(3.6) 


and 


0  +  Aa  = 


Subtract  Eq.  (3.6)  from  (3.7)  to  obtain 


tn 


/r  +  Ar\ 
'  ro  / 


(3.7) 


Aa 


=  - - i -  *n(l  +  ~) 

c^/g(f')max  V  ' 


Solving  for  Ar  we  find  that 


where 


sr  =  r0(eQ-l) 


(3.8) 
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We  can  generalize  Eq.  (3.8)  for  any  step  size  Arn  by  replacing  rQ  with 

rn _i ,  the  last  propagation  distance.  On  the  basis  of  Pestorius's  expe- 
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rience  we  choose  Aa=0.1.  We  use  a  subroutine  NWSLOPE  to  calculate  the 
slope  of  the  head  shock  (the  steepest  part  of  the  waveform).  The  head 
shock  of  the  input  waveform  usually  contains  only  three  or  four  points. 
Therefore,  we  use  a  scheme  in  which  we  find  the  slope  over  the  central 
portion  of  the  head  shock.  We  first  decide  between  which  two  points  we 
want  to  compute  the  slope.  The  points  are  located  at  certain  fractions  of 
the  N  wave  amplitude  (AMAX).  The  fractions  (FI  and  F2)  are  input  from  the 
terminal.  The  time  increments  (T1  and  T2)  corresponding  to  (Fl-AMAX) 
and  F2-AMAX)  are  found  by  performing  a  linear  interpolation  between  the 
two  points  bounding  (FI -AMAX)  and  (F2-AMAX),  respectively.  The  slope  is 
then  calculated  using  the  following  equation: 

SLOPE  =  ( F2-F1 )  •  AMAX/(T2-T1 )  .  (3.9) 

We  can  now  write  Eq.  (3.8)  in  code.  The  variable  AXSTP,  in  code, 

2  / 

corresponds  with  Aac/6.  Q  now  becomes 

Q  *  AXSTP/ (RN- SLOPE) 

Finally  the  step  size  increment  (SSI)  is 

SSI  =  RN(eQ-l)  .  (3.10) 

c.  Application  of  Absorption 

We  transform  from  the  time  domain  to  the  frequency  domain  to 
apply  absorption  and  dispersion  to  the  waveform.  The  absorption  we 
consider  is  actual  atmospheric  absorption  (quiet  air).  The  Applied 
Research  Laboratories  computer  library  subroutine  FFTCC  was  used  to 
transform  to  the  frequency  domain.  The  absorption  and  dispersion  are 
calculated  by  the  function  DABSORP,  which  is  a  modified  version  of  the 
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function  ABSORP,  also  to  be  found  in  the  user  library.  It  is  a  digital 
implementation  of  ANSI  document  SI. 26-1978,  "Standard  Method  for  the 
Calculation  of  the  Absorption  of  Sound  by  the  Atmosphere". ^  This  function 
allows  one  to  compute  the  absorption  coefficient  for  a  given  frequency, 
ambient  temperature,  ambient  pressure,  and  relative  humidity.  A  more 
detailed  discussion  of  the  function  is  found  in  the  discussion  of  the 
subroutines. 

The  absorption  coefficients  are  calculated  in  the  following 
manner.  The  number  of  frequencies  at  which  we  compute  absorption 
coefficients  is  NF=N/2,  where  N  is  the  number  of  points  in  the  waveform. 
There  is  no  absorption  coefficient  for  the  dc  component.  We  compute  only 
N/2  coefficients  because  the  highest  frequency  we  can  discern  is  half  the 
sampling  frequency  (fs=l/DELT).  The  frequency  increment  is 

DELF  =  N  *VlT 

where  DELT  is  the  time  increment  between  points.  The  frequencies  at  which 
the  absorption  coefficients  are  calculated  are  found  by  multiplying  DELF 
times  J,  where  J  is  a  counter  from  1  to  NF.  The  value  for  the  absorp¬ 
tion  coefficient  (in  dB/m)  at  each  frequency  is  then  computed  by 
DABSORP.  Thermovi scous  absorption,  rotational  relaxation  of  nitrogen  and 
oxygen,  and  vibration  re’axation  of  nitrogen  and  oxygen  are  included  in  the 
absorption  coefficient.  In  addition  to  the  absorption  coefficient,  the 
function  returns  the  value  of  the  oxygen  relaxation  frequency  and  the 
absorption  coefficient  due  to  oxygen  relaxation.  We  find  the  quantity  n, 
which  is  related  to  the  dispersion  (see  Eqs.  (3.1a)  and  (3.1b)),  from 
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these  values.  Since  the  absorption  coefficient  is  in  dB/m  we  need  to 
modify  Eq.  (3.1 ) , 


U  = 


r 

o 


r  +  Ar 
o 


-c(r-r  )/ 20 
10  0 


(3.11) 


For  our  case 


a  =  aTV  +  aN2 


) 


In  the  program  the  complex  variable  (-0  is  named  ALFA  and  the  spherical 
spreading  factor  ^/(r^+Ar)  is  named  SPREAD.  Eq.  (3.11)  is  written  in 
code  as 

A(J)  =  SPREAD/N  *  A(J)old  i0(ALFA*SSI/20) 

The  factor  1/N  is  needed  for  normalization  because  of  certain  properties 
of  the  subroutine  FFTCC. 

d.  Application  of  Nonlinear  Distortion 
After  the  absorption  step  is  taken  we  transform  back  to  the 
time  domain  to  apply  the  nonlinear  distortion.  The  distortion  is  applied 
using  Eq.  (3.5).  In  our  program  the  term  BAr/CQ  is  named  DISFAC. 

Equation  (3.5)  is  represented  in  code  as 

T(J)  =  C(J-1 )  *  DELT]  -  [DISFAC  *  P(J)] 

old- 


where  ( (J-l )*DELT)=T(J) 
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After  the  distortion  step,  the  time  increment  between  points  is  no  longer 
the  same  throughout  the  waveform.  The  FFT  routine  we  use  requires  data 
sampled  at  equal  time  increments.  To  resample  the  data  points  at  equal 
time  increments  we  used  a  subroutine  called  RESAMP.  This  subroutine  is 
discussed  in  the  section  describing  the  subroutines, 
e.  Data  Output 

The  program  returns  to  the  step  size  calculation  after  re¬ 
sampling  the  data  points.  When  the  propagation  distance  reaches  one  of 
the  specified  propagation  distances,  the  time  waveform  is  plotted  by  means 
of  the  subroutine  NWPLOT.  If  a  frequency  spectrum  has  also  been  speci¬ 
fied,  it  is  plotted  at  this  time  using  the  subroutine  FQPLOT.  A  plot  file 
is  created  in  both  plot  subroutines.  These  subroutines  are  discussed  in 
the  next  section. 

The  entire  process  is  repeated  until  the  waveform  reaches  its 
final  propagation  distance.  At  the  end  of  the  program  the  plot  file  is 
sent  to  a  Nicolet  Zeta  Research  plotter,  which  draws  the  waveforms. 

2.  Subroutines 
a.  DABSORP 

The  function  DABSORP  is  a  modified  version  of  the  Applied 
Research  Laboratories  computer  library  function  ABSORP.  DABSORP  is 
accurate  to  within  ±10%  over  the  temperature  range  255. 4°K  (0°F)  to 
310. 9°K  (100°F)  at  relative  humidities  from  0%  to  100%,  frequencies 
100  Hz  to  10  MHz,  and  ambient  pressures  from  0  to  2  atm.  It  was 
necessary  to  modify  ABSORP  to  include  the  effect  of  dispersion  due  to 
oxygen  relaxation.  Our  modification  utilizes  a  separate  calculation  of 
the  absorption  coefficient  due  to  oxygen  relaxation  (0ALPH).  The 
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modified  function  passes  the  value  of  OALPH  as  well  as  the  value  for  the 
relaxation  frequency  of  oxygen  (FRO)  to  the  driver  program.  The  absorp¬ 
tion  coefficient  used  in  the  absorption  step  (ALFA)  is  defined  as 

ALFA  =  -e 

=  -a  +  jn 

=  ALR  +  j(ALI)  , 

where 

ALR  =  -a  =  -  DASORP,  and 

ALI  =  n  =  ioia  =  a 
0 

=  F  *  OALPH/FRO 

F  is  the  frequency  at  which  we  calculate  the  absorption  coefficient.  A 
listing  of  DABSORP  is  found  in  Appendix  B. 

b.  RESAMP 

The  nonlinear  distortion  step  shifts  the  time  values  of  the  data 
points  by  an  amount  proportional  to  their  particle  velocity.  The  sub¬ 
routine  RESAMP  uses  a  point- si  ope  formula  to  interpolate  between  the 
unequally  spaced  time  values,  to  find  values  at  equal  time  intervals,  so 
that  the  FFT  routine  can  be  used.  The  subroutine  call  is 
CALL  RESAMP  (P,T,WK,NR,N,DELT,IERR) . 

WK  is  a  scratch  array  of  dimension  NR.  IERR  is  an  error  flag  which 
returns  1  if  a  sampling  error  is  encountered  and  0  on  no  error.  A  list¬ 
ing  of  RESAMP  is  found  in  Appendix  B. 

c.  NWPLOT 


This  subroutine  is  used  to  plot  time  waveforms  at  the  desired 
propagation  distances.  For  each  run,  the  amplitude  of  each  waveform  is 
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normalized  to  the  amplitude  of  the  input  waveform.  If  the  ratio 
P/PNORM  is  small,  it  is  multiplied  by  an  appropriate  scale  factor.  As  a 
result,  all  of  the  waveforms  (even  those  at  large  propagation  distances) 
are  large  enough  for  us  to  analyze  accurately.  The  subroutine  call  is 
CALL  NWPLOT  (P,T,N,PNORM,IRESULT,TMIC,RN) 

PNORM  is  the  amplitude  of  the  input  waveform.  IRESULT  is  the  result 
file  which  returns  a  listing  of  the  P  and  T  arrays  at  each  specified 
propagation  distance.  TMIC  is  the  number  of  vsec/in.  specified  for 
X  axis.  RN  is  the  propagation  distance.  A  listing  of  NWPLOT  is  found 
in  Appendix  B. 

d.  EQPLOT 

FQPLOT  plots  the  frequency  spectrum  of  the  N  wave.  The 
spectrum  is  plotted  on  a  log-log  scale  normalized  to  the  peak  particle 
velocity  so  that  at  the  fundamental  frequency  the  amplitude  is  0  dB. 

Use  of  this  subroutine  is  optional.  Selection  is  made  in  the  driver 
program  by  picking  IOPT  equal  to  0  (frequency  plot  only)  or  2  (frequency 
plot  and  time  waveform).  The  subroutine  call  is 

CALL  FQPLOT  (A,DELT,N,IFIRST) 

I  FIRST  is  a  flag  that  signals  the  first  propagation  step.  A  listing  of 
FQPLOT  can  be  found  in  Appendix  B. 

3.  Frequency  Domain  Analysis 

In  an  effort  to  determine  the  effect  of  dispersion  of  the  N 
wave,  we  performed  many  computer  propagations,  both  with  and  without 
dispersion.  We  found  very  little  difference  between  the  time  waveforms 
that  were  computed  with  dispersion  and  those  that  were  computed  without 
diversion.  The  effect  of  dispersion  on  a  finite  amplitude  sine  wave 
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has  been  studied  by  Blackstock.  In  order  to  check  whether  the 
dispersion  was  working  correctly  in  our  program,  we  input  a  sine  wave 
and  computed  the  amplitude  of  the  fundamental,  second  harmonic,  and  third 
harmonic.  Our  results  did  not  agree  with  those  predicted  by  Blackstock. 
We  discovered  that  the  FFT  routine  was  the  cause  of  our  discrepancy.  As 
noted  in  the  section  describing  the  application  of  absorption,  the  number 
of  filters  or  frequency  components  depends  on  the  number  of  points  in  the 
waveform.  These  filters  have  a  (sinx)/x  weighting  function  associated 
with  them.  Therefore,  the  amplitudes  returned  by  the  program  are  not 
for  discrete  frequencies  but  for  a  frequency  "band".  To  get  a  more 
accurate  measurement  of  the  amplitudes  we  employed  an  aperture  shading 
subroutine  called  APWT  in  which  Taylor"^  weighting  was  used.  The  peak- 
to-si delobe  ratio  was  30  dB.  After  the  routine  was  incorporated  in  the 
program,  we  found  that  the  computed  amplitudes  agreed  well  with  the 
analytical  predictions.  As  a  further  check  we  ran  the  program  with  and 
without  the  weighting  to  determine  whether  the  time  waveforms  were 
affected  as  much  as  the  frequency  spectrum.  We  found  the  difference 
between  the  runs  was  negligible,  and  concluded  that  whatever  was  done  in 
the  FFT  was  undone  by  the  (FFT)""'.  It  was  concluded  that  the  aperture 
shading  was  needed  only  if  one  desires  a  frequency  analysis  of  the  wave¬ 
form.  Since  only  time  waveforms  were  analyzed,  we  decided  to  delete  the 
aperture  weighting  subroutine.  However,  a  listing  of  APWT  and  TAYLOR  is 
found  in  Appendix  B,  and  should  be  included  if  a  frequency  analysis  is 
done. 
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D.  Results 

We  obtained  very  good  agreement  between  our  experimental  wave¬ 
forms  and  the  computed  generated  waveforms.  The  results  are  discussed 
in  more  detail  in  Chapter  5.  Two  sample  output  waveforms  are  shown  in 
Fig.  3.2.  The  amplitudes  of  the  waveforms  are  found  in  the  result  file 
generated  by  the  driver  program.  The  value  m^  is  the  maximum  slope  of 
head  shock.  The  maximum  slope  for  our  measured  waveforms  is  consistently 
the  slope  of  the  tangent  to  the  head  shock  at  p  =  0.5  P.  Therefore,  we 
adopt  this  method  for  measuring  our  waveforms.  The  rise  time  is  simply 
the  amplitude  P  divided  by  m-j .  The  half  duration  T  is  the  time  between 
the  foot  shock  and  the  zero  crossing  of  the  expansion  phase.  We  have 
defined  the  shock  foot  to  occur  at  the  intersection  of  the  tangent  line 
with  the  time  axis.  This  definition  is  consistent  with  our  experimental 


observations. 


CHAPTER  4 


EXPERIMENT 

This  chapter  contains  a  detailed  description  of  the  measurement 
apparatus  and  method.  A  few  sample  measurements  are  shown. 

A.  Apparatus 

The  measurements  system  has  four  main  parts:  spark  source, 
optical  bench,  receiving  system,  and  data  capture  and  storage.  The 
arrangement  of  the  system  is  shown  in  Fig.  4.1.  The  important  features  of 
each  are  described  in  the  following  paragraphs. 

1 .  Spark  Source 

A  Grozier  Technical  Systems  Incorporated  Model  51-4  spark 
generator  was  used  to  provide  energy  for  the  spark  by  charging  a  1.0  pF 
capacitor  to  0-4  kv.  The  capacitance  was  allowed  to  discharge  across  a 
gap  between  two  diametrically  opposed  tungsten  electrodes.  Typical  gap 
lengths  were  between  1  and  5  mm. 

A  third  electrode,  located  perpendicular  to  the  gap,  was  used 
to  trigger  the  discharge  at  a  set  voltage.  The  advantages  of  this  system 
are  flexibility  and  reproducibility.  Without  the  trigger  electrode  we 
are  restricted  to  one  voltage,  the  natural  breakdown  potential,  for  each 
gap  length.  Moreover,  the  natural  breakdown  potential  varies  from  spark 
to  spark.  With  the  trigger  in  use,  when  the  voltage  across  the  capacitor 
reaches  the  voltage  preset  by  the  operator,  a  small  spark  jumps  from  the 
3rd  electrode  and  stimulates  breakdown.  The  operator  is  therefore  able 
to  control  the  spark  energy,  which  largely  determines  amplitude  and  half 


45 


46 


MICROPHONE 

AMPLIFIER  PREAMPLIFIER 


SPARK 


LASER 


<  t 


MICROPHONE 


LASER  BEAM 
AND  AXIS  OF 
SYMMETRY 


TRIGGER 
ANTENNA i 


OPTICAL  BENCH 


PULSE  GENERATOR 


NICOLET 

OSCILLOSCOPE 


TRIGGER 


X-Y  RECORDER 


FIGURE  4.1 

EXPERIMENTAL  ARRANGEMENT 


•'  » 

t  ■  i 


OPTICAL  BENCH 
CROSS  SECTION 


ARL:UT 

AS-82-122 

LBO-GA 

2-1682 


47 

duration  of  the  N  wave.  The  breakdown  voltage  was  monitored  by  a  Simpson 


Model  260  voltmeter. 

The  electrode  configuration  is  shown  in  Fig.  4.2.  Each 
electrode  was  1/16  in.  in  diameter  and  its  point  was  ground  to  a  cone  of 
half-angle  approximately  30°.  A  polypropylene  block  was  machined  to 
support  the  electrodes.  The  electrodes  were  held  in  threaded  aluminum 
dowels.  The  spark  gap  could  be  adjusted  symmetrically  by  turning  the 
dowels.  A  photograph  of  an  actual  spark  discharge  is  found  in  Fig.  4.3; 
neutral  density  filters  were  used  to  cut  down  on  the  light  generated  by  the 
spark.  It  is  clear  that  the  breakdown  does  not  really  take  place  along  a 
specific  path  but  in  a  general  region  around  the  electrodes. 

2.  Optical  Bench 

The  I-beam  rails  used  in  this  experiment  were  designed  by 
30 

Anderson  and  are  described  in  his  thesis.  They  were  made  in  the  machine 
shop  at  Applied  Research  Laboratories.  Two  rails  were  used  as  an  optical 
bench;  the  first  one  12  ft  long  and  the  second  one  8  ft  long.  The  rails 
were  clamped  on  a  heavy  table  and  leveled  to  within  0.0005  in. /ft.  A 
Hughes  Model  3193H  laser  was  used  to  align  the  tracks  of  the  two  rails. 

3.  Receiving  System 

The  condenser  microphone  used  was  based  on  designs  and 
construction  techniques  developed  over  the  years  by  previous  students  at 

Applied  Research  Laboratories.  The  first  microphones  were  built  by 

35  27  30 

Cornet,  who  modified  designs  by  Wright.  Anderson  improved  the 

construction  process  and  was  able  to  obtain  a  high  degree  of  consistency 
from  one  microphone  to  the  next.  Basically,  a  condenser  microphone  con¬ 
sists  of  a  metallic  diaphragm,  which  is  exposed  to  the  sound  field,  a 


FIGURE  4.3 
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compressible  dielectric  layer  (commonly  air),  and  a  metallic  backplate. 

The  diaphragm  and  the  backplate,  separated  by  the  dielectric,  act  as  a 
condenser,  across  which  is  imposed  a  dc  bias  voltage  (typically  about 
200  V).  Sound  pressure  on  the  diaphragm  changes  the  spacing  between  it 
and  the  backplate,  thus  altering  the  capacitance  of  the  microphone.  If  a 
constant  charge  is  maintained  on  the  condenser  plates  a  corresponding 
voltage  is  induced,  which  is  fed  to  the  input  of  a  preamplifier. 

Cornet's  and  Anderson's  microphones  were  built  utilizing  the 
case  of  a  standard  BNC  connector  as  an  outer  housing.  A  plexiglass 
insulator  was  machined  and  force- fit  or  cemented  into  the  outer  housing. 

A  small  hole  was  drilled  in  the  insulator  for  the  backplate  electrode, 
which  was  then  fit  into  the  hole  and  cemented  in  place.  The  diaphragm 
material  used  was  1/8  mil  mylar  metal i zed  on  one  side.  The  mylar  was 
laid,  metalized  side  up,  directly  on  the  backplate  surface.  The  only  air 
gap  (and  only  stiffness)  was  provided  by  tiny  pockets  of  air  trapped  in 
the  microscopically  rough  surface  of  the  backplate. 

Diffraction  at  the  outer  edges  of  the  microphone  face  gives  rise 
to  an  inverted  replica  of  the  incident  wave.  If  the  outer  housing  is 
small,  the  diffracted  signal  arrives  shortly  after  the  incident  wave  and 
interferes  with  it.  If  the  microphone  housing  is  made  sufficiently  large, 
however,  the  diffracted  signal  is  delayed  long  enough  that  it  may  be 
ignored.  Cornet  and  Anderson  loosely  fit  a  plane,  plexiglass  baffle  about 
the  microphone  housing  to  delay  the  diffracted  signal.  Our  first  rise 
time  measurements  were  made  with  this  system.  However,  it  was  very 
difficult  to  fit  the  baffle  securely  on  the  microphone  housing.  There 
tended  to  be  a  gap  or  step  at  the  baffle-microphone  juncture.  The 
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discontinuity  at  the  juncture  apparently  gave  rise  to  diffracted  waves, 
which  contaminated  the  received  signal. 

We  modified  the  microphone  design  hoping  to  avoid  signal 
contamination.  New  microphones  combining  the  microphone  and  baffle  in 
one  system  were  built.  The  diameter  of  the  microphone  housing  was 
expanded  to  1.5  in.  The  remainder  of  the  front  surface  of  the  microphone 
is  extra  insulator  area  to  act  as  a  baffle.  A  cutaway  view  of  our  micro¬ 
phone  is  shown  in  Fig.  4.4(a).  The  outer  housing  is  machined  from  brass 
or  aluminum.  A  BNC  connector  is  attached  beneath  the  housing.  Cornet's 
method  of  fitting  the  insulator  and  backplate  in  the  housing  was  used. 

The  procedure  used  to  finish  the  backplate  is  the  most  critical 
step  in  the  construction  of  the  microphone.  First,  a  lathe  was  used  to 
machine  the  face  of  the  microphone  flat.  Next,  the  backplate  was  smoothed 
using  No.  600  silicon  carbon  sandpaper.  The  microphone  was  then  placed 
in  an  ultrasonic  cleaner  to  remove  the  grit  particles.  Finally,  the 
surface  was  buffed  and  polished.  This  process  created  a  random  distribu¬ 
tion  of  air  filled  microscopic  cavities,  each  having  a  different  resonance 
frequency.  The  mass  of  the  diaphragm  of  each  cavity  and  the  spring  formed 
by  the  air  trapped  in  the  cavity  is  modeled  as  a  small  oscillator.  The 
resonance  frequency  of  the  microphone  is  determined  by  the  overall 
resonance  of  all  the  cavities  operating  in  parallel.  The  response  of  many 

oscillators  in  parallel,  each  having  a  different  resonance  frequency,  is 
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much  flatter  than  that  for  a  single  oscillator;  see  Cornet.  As  a 
result,  the  microphone  is  less  likely  to  ring.  To  make  the  resonance 
frequency  of  each  individual  oscillator  as  high  as  possible,  the  oscil¬ 
lator  mass  should  be  as  small  as  possible,  i.e.,  the  mylar  should  be  as 


ALUMINIZED 
MYLAR  (0.1  mil) 


52 


SILVER 


(a)  CUTAWAY  VIEW  OF  THE  MICROPHONE 


(b)  FIELD  EFFECT  TRANSISTOR  PREAMPLIFIER 


FIGURE  4.4 

MICROPHONE  AND  PREAMPLIFIER 


ARLUT 
AS- 82-1 28 
LBO  -  Q  A 
2- 16-82 


l 

« 


53 

thin  as  possible.  The  mylar  used  for  our  microphones  was  1/10  mil 
(Sheldahl  Company,  Northfield,  Minnesota  55057). 

We  designed  the  backplate  such  that  the  measurement  of  rise 
time  would  not  be  limited  by  the  spherical  waves  we  intend  to  use.  A 
backplate  diameter  of  1/8  in.  was  chosen.  The  difference  between  the 
time  a  spherical  wave  of  radius  10  cm  reaches  the  center  of  the  back¬ 
plate  and  the  time  it  reaches  the  edge  of  the  backplate  is  0.04  usee.  We 
assume  that  the  microphone  is  effectively  a  point  detector  because  the 
active  area  is  very  small.  In  particular,  if  the  source  and  microphone 
are  misaligned  by  1°,  the  response  to  a  1.6  MHz  signal  is  down  0.4  dB. 
Below  56  kHz  the  microphone  is  omnidirectional  (eHp=90°)  within  ±3  dB. 

The  FET  preamplifier  constructed  by  Cornet  was  used  in  our  early 
experiments.  We  found  that  the  dynamic  range  was  insufficient  for  measur¬ 
ing  strong  signals  close  to  the  spark.  L.  A.  Thompson  designed  a  new 
preamplifier,  similar  in  design  to  Cornet's,  that  has  a  much  greater 
dynamic  range  and  a  better  frequency  response.  A  circuit  diagram  is  shown 
in  Fig.  4.4(b).  The  maximum  voltage  swing  is  10  V  peak  to  peak.  The  low 
and  high  frequency  3  dB  down  points  are  10  Hz  and  2.25  MHz,  respectively. 
The  input  capacitance  is  9.9  pF,  approximately  one-fifth  that  of  Cornet's 
preamplifier.  The  preamplifier  contains  the  bias  voltage  (180  V)  for 
the  microphone.  Also  included  is  a  line  driver,  which  drives  a  50  n 
cable  without  much  loss.  The  switch  shown  is  used  to  disconnect  the 
batteries  so  they  do  not  drain  while  the  preamplifier  is  not  in  use.  An 
early  difficulty  with  the  preamplifier  was  a  dc  offset  in  the  output.  A 
100  kn  potentiometer  was  installed  in  the  circuit  to  allow  the  offset  to 
be  balanced  out. 
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In  increasing  the  frequency  response  of  the  microphone- 
preamplifier  system,  we  lost  quite  a  bit  of  sensitivity.  To  avoid  the 
cable  noise  interfering  with  the  signal,  it  was  necessary  to  amplify  the 
signal  before  sending  it  to  the  oscilloscope.  We  used  an  HP461A 
amplifier  for  this  purpose. 

There  were  some  limitations  on  our  receiving  system.  The 
HP461A  amplifier  has  a  low  frequency  3  dB  down  point  of  360  Hz.  The 
FET  preamplifier  has  a  high  frequency  3  dB  down  point  of  2.25  MHz, 
and  the  microphone  has  a  high  frequency  3  dB  down  point  of  about  1  MHz. 
The  rolloff  at  the  low  frequency  end  is  determined  by  the  HP461A  ampli¬ 
fier  (18  dB/octave).  The  FET  preamplifier  rolls  off  at  the  high  fre¬ 
quency  end  at  6  dB/octave.  The  high  frequency  rolloff  of  the  system, 
however,  is  set  by  the  condenser  microphone.  Generally,  a  condenser 
microphone  rolls  off  at  12  dB/octave.  However,  since  our  microphone  is 
really  many  small  microphones,  each  with  a  different  resonance  frequency, 

in  parallel,  the  rolloff  is  expected  to  be  much  gentler. 
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Crocker  has  developed  a  model  for  the  effect  of  limited  band¬ 
width  of  a  receiver  on  response  to  an  ideal  N  wave.  The  model  is  valid 
only  for  rolloffs  of  6  dB/octave.  Therefore  we  cannot  adjust  our  wave¬ 
forms  using  the  results.  However,  we  use  his  model  to  get  a  general  idea 
of  the  effects  of  our  limited  bandwidth.  Crocker  dealt  with  the  low 
frequency  limitation,  which  can  be  modeled  as  a  high  pass  filter,  and  the 
high  frequency  limitation,  which  can  be  modeled  as  a  low  pass  filter 
separately.  His  analysis  is  shown  in  Fig.  4.5.  In  the  case  of  low 
frequency  limitation,  the  effects  vary  depending  on  the  ratio  a=T-|/T0, 
where  T-|  is  the  RC  time  constant  of  the  high  pass  filter.  As  a  decreases 
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From  the  raport  by  M. J.  Crocker,  "Maaauramant  of  Sonic  Boom  with  Limited  Frequency  Retponae 
I netru mentation:  A  Theoretical  Study,"  Wyle  Labi  Rawarch  Staff  Report  WR  66-20,  April  1966. 


from  infinity  to  zero,  the  midsection  of  the  wave  distorts,  thus  decreas¬ 
ing  the  half  duration  (see  Fig.  4.5(a)).  Using  the  model  to  estimate  the 
reduction  in  half  duration  for  our  system  we  found  that  the  maximum  error 
is  4.1%.  The  high  frequency  limitation  is  determined  by  the  ratio  6=T2/Tq 
where  is  the  RC  time  constant  of  the  low  pass  filter.  As  6  grows 
from  zero  to  infinity,  the  rise  time  increases,  amplitude  decreases,  and 
the  half  duration  increases  (see  Fig.  4.5(b)).  For  our  high  frequency 
3  dB  down  point  (1  MHz),  we  find  that  T~=l  usee.  For  an  N  wave  with 
To=20  psec  (a  typical  value  for  our  experimental  N  waves),  this  gives 
6=1/20=0.05.  By  interpolating  between  the  6=0.02  and  6=0.1  curves,  one 
can  see  that  a  rise  time  of  order  tr  =  0.2  To  =  4  psec  is  indicated. 
However,  we  were  able  to  measure  rise  times  as  low  as  1  psec.  Therefore, 
we  conclude  that  Crocker's  model  for  the  high  frequency  limitation  does 
not  fit  our  receiving  system. 

4.  Data  Capture  and  Storage 

The  signal  from  the  preamplifier,  after  going  through  the  HP461A 
amplifier,  was  input  to  a  Nicolet  Model  2090- III  digital  oscilloscope 
equipped  with  a  Model  206-2  plug-in  unit.  A  magnetic  disk  memory  was  used 
to  store  signals  for  further  analysis.  The  disk  is  divided  into 
8  "tracks",  each  having  a  capacity  of  4096  data  points.  Each  "track"  can 
be  divided  into  "quarters",  each  having  a  capacity  of  1024  data  points. 

Use  of  an  entire  "track"  permitted  the  storage  of  2  msec  of  data  with  a 
time  resolution  of  0.5  psec/point  and  12-bit  accuracy.  Electromagnetic 
radiation  was  picked  up  by  a  loop  antenna,  amplified  by  a  HP450A  amplifier 
and  then  used  to  open  the  gate  of  a  Chronetics,  Inc.,  Model  PG-11  pulse 
generator,  which  in  turn  triggered  the  oscilloscope.  A  delay  set  by  the 


width  of  the  pulse  from  the  Chronetics  unit  made  it  possible  to  hold  the 
oscilloscope  trigger  until  the  acoustic  wave  arrived  at  the  microphone. 
Captured  waveforms  were  either  saved  on  the  disk  or  plotted  directly  on 
HP7Q45A  X-Y  recorder. 

B.  Calibration 

The  technique  used  to  calibrate  the  microphone-preamplifier 

37 

system  has  been  described  by  Davy  and  Blackstock  and,  in  more  detail, 

35  30  38  32 

by  Cornet,  Anderson,  Cobb,  and  Essert.  The  technique  is  based 

upon  the  lengthening  and  excess  attenuation  of  an  outgoing  spherical  N 

wave.  We  use  only  the  positive  half  of  our  spark  produced  N  waves, 

because  the  negative  half  is  not  sufficiently  close  to  the  ideal  N  wave 

shape.  The  peak  pressure  amplitude  P  and  half  duration  T  of  an  ideal, 

spherically  spreading  N  wave  at  radial  distance  r  are  related  to  the 

pressure  PQ  and  half  duration  Tq  at  some  reference  position  rQ  by  the 

following: 

*  -  ■'oPo(1  +  °o  tn(y]  -  <4-'l 

T  ■  To['  +  "o  en(yV2  •  l4-2) 

where 

=  (r±L)„  J>P£  (4  3) 

2p  c  o 
o  o 


Equation  (4.1)  is  multiplied  by  Eq.  (4.2)  to  find 

TrP  =  VoPo  • 


(4.4) 


which  shows  that  the  triple  product  has  the  same  value  regardless  of  the 
distance.  The  peak  pressure  P  and  the  peak  voltage  E  are  related  by 


P  =  E/S 

where  S  is  the  sensitivity  of  the  microphone-preamplifier  system. 
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(4.5) 


Combining  Eqs.  (4.3),  (4.4),  and  (4.5)  yields 


S  = 


O  O  O  0 

Ty+1T 


(4.6) 


We  have  replaced  TQroEo  by  TrE,  an  average  value  of  the  triple  product, 
in  Eq.  (4.6)  because  the  accuracy  of  an  average  is  greater  than  that  of 
a  single  measurement.  The  parameter  TrE  is  determined  from  measurements 
of  the  N  wave  amplitude  and  half  duration  at  various  distances.  If 
Eq.  (4.2)  is  squared,  the  result  is 


T2  =  T2 


n  +  °nTn  WJM 
O  0  0  \r0 ) 


2  2 
The  quantity  aQTQ  may  thus  be  found  by  plotting  data  for  T  against 

tn(r/r  )  and  measuring  the  slope  of  the  straight  line  that  results.  A 

least  squares  fit  is  used. 

Although  with  the  Grozier  spark  generator  a  spark  could  be 

obtained  at  a  consistent  breakdown  voltage,  it  was  found  that  T  and  P 

varied  slightly  with  each  spark.  To  reduce  this  variability  we  used  a 

32 

special  procedure  developed  by  Essert.  An  "average"  waveform  was  found 
as  follows:  At  each  measurement  position  ten  N  waveforms  were  captured 
and  stored  on  the  disk.  The  values  for  T  and  P  were  recorded  for  each 
waveform.  Average  values  for  T  and  P  were  found  from  the  data.  Then  we 
captured  several  waveforms  until  we  found  one  that  had  the  average  values 

of  T  and  P.  This  waveform  was  plotted  on  the  X-Y  recorder  and  used  for 

further  analysis. 
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Ringing  and  a  finite  rise  time  were  found  in  the  measured 

waveforms.  It  was  necessary  to  idealize  the  measured  waveforms,  since 

the  calibration  procedure  is  based  on  an  ideal  N  wave.  The  ringing  in 

the  linear  decay  portion  of  the  waveform  was  smoothed  out  by  assuming 

that  the  actual  waveform  passed  through  the  center  of  oscillations.  We 

placed  the  idealized  head  shock  at  a  point  midway  between  the  foot  and 

the  peak  of  the  measured  waveform.  The  procedure  is  shown  in  Fig.  4.6. 

The  peak  voltage  E  and  half  duration  T  were  determined  from  the  ideal 

waveform  and  used  in  the  calculations  for  calibration.  For  further 
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descriptions  on  the  idealizing  procedure,  see  Cobb,  Anderson,  and 
32 

Essert. 

Two  different  microphone-preamplifier  systems  were  used  in  our 

experiment.  The  first  system  consisted  of  one  of  the  small  microphones, 

with  fitted  baffle,  and  Cornet's  preamplifier.  N  waves  were  recorded  at 

eight  source- receiver  distances  from  10  to  180  cm  for  a  spark  energy  of 

1.04  J.  The  temperature  was  75° F  and  the  relative  humidity  66%. 

_  2 

Figure  4.7  shows  the  calculation  of  TrE  and  aoTo.  The  results  are 
summarized  as  follows 

TrE  =  144.4  ysec  x  cm  x  v,  standard  deviation  =  3.26%, 

2  2 

oqT0  =  20.223  nsec  ,  correlation  coefficient  =  0.996, 
and,  finally, 

S  =  5.68  mbar/V 


or 


S  =  -75.09  dB  re  1  V/ybar, 


! 


FIGURE  4.6 

IDEAL  N  WAVE  (ABOVE)  AND  IDEALIZATIONS  OF  MEASURED  N  WAVES  (BELOW) 
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where  S  =  20  log  and  Sref  is  the  reference  sensitivity  1  V/ubar. 

This  value  is  the  sensitivity  of  the  entire  microphone-preamplifier 
system. 

The  second  system  consisted  of  the  new  preamplifier  and  the  new, 

larger  microphone.  In  the  calibration  run  for  this  system  N  waves  were 

recorded  at  ten  source-receiver  distances  from  30  to  200  cm  for  a  spark 

energy  of  4.35  J.  The  temperature  was  73° F  and  the  relative  humidity  43%. 

_  2 

The  methods  described  above  were  used  to  find  TrE  and  a  T  .  The  results 

o  o 

are  summarized  as  follows  (see  Fig.  4.8). 

TrE  =  27.357  ysec  *  cm  x  V,  standard  deviation  =  6.9% 

2  2 

a  T  =  151.92  usee  ,  correlation  coefficient  =  0.987, 

0  o 

and,  finally 

S  =  225.23  mbar/V 
or 

S  =  -107.05  dB  re  1  V/ubar. 

This  is  the  value  of  the  sensitivity  of  the  microphone-preamplifier 
system,  with  0  dB  gain  from  the  HP461A  amplifier. 

Since  the  calibration  method  used  was  an  unconventional  one, 
a  completely  different  calibration  was  performed  to  check  our  results. 

We  used  a  comparison  method,  in  which  the  same  sound  field  is  measured 
with  two  different  microphones,  one  of  which  has  a  known  calibration. 

The  other  microphone  was  our  new  one.  A  plane  wave  field  in  a  2  in.  diam 
tube  was  chosen  to  make  this  calibration.  The  tube  was  3.283  ft 
long.  A  JBL  Model  375  horn  driver  was  attached  on  one  end  of  the  tube. 

On  the  other  end  we  attached  a  flange,  specially  machined  to  hold  our 
microphone  (designated  microphone  MM)  and  a  Bruel  and  Kjaer  (B&K)  1/4  in. 


’  »  » 
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microphone.  A  cutaway  view  of  the  flar  a  is  shown  in  Fig.  4.9.  Both 
microphones  were  mounted  flush  with  the  inner  wall  of  the  flange.  Pulsed 
sine  waves  at  test  frequencies  of  0.5,  1,  2,  and  4  kHz  were  sent  down  the 
tube  via  the  driver,  and  the  signal  was  picked  up  with  both  microphones. 
Prior  to  the  calibration  run  the  B&K  microphone  was  calibrated  with  a 
pistonphone.  The  output  signals  from  the  two  microphones  were  captured 
on  the  Nicolet  oscilloscope,  which  was  triggered  by  the  pulser,  and 
their  amplitudes  recorded  in  millivolts  (see  Fig.  4.10  for  setup).  The 
sensitivity  of  microphone  MM  was  calculated  as  follows. 


=  20  log 


(4.7) 


where  SgK  is  the  sensitivity  of  the  B&K  microphone,  Sref  is  the  reference 
sensitivity  (1  V/ybar),  EMM  is  the  output  voltage  of  microphone  MM,  and 
EgK  is  the  output  voltage  of  the  B&K  microphone.  A  plot  of  sensitivity 
versus  frequency  gives  the  frequency  response  of  the  microphone- 
preamplifier  system  (see  Fig.  4.11).  Since  we  expected  the  frequency 
response  to  be  flat  in  the  frequency  region  that  we  measured  (0.5-4  kHz), 
the  best  measure  of  sensitivity  was  an  average  of  the  values  at  the  four 
test  frequencies. 

Since  the  two  microphones  were  physically  separated,  they 
measured  slightly  different  sound  fields.  We  account  for  this  by  adjusting 
the  amplitude  recorded  by  the  B&K  microphone  to  the  value  it  would  have 
been  if  it  had  been  in  exactly  the  same  position  as  microphone  MM.  If 
the  incident  wave  at  microphone  MM  is  described  by 
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FLANGE  FOR  CALIBRATION  IN  2  in.  TUBE 
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I  =  A  eju>t 

then,  because  the  surface  is  rigid,  the  reflected  wave  also  has  an 
amplitude  A  so  that  the  total  pressure  is 

P  =  2A 
kMM  e 

The  combination  of  the  incident  plus  reflected  signals  seen  by  the  B&K 

microphone,  which  is  Ax  away  from  microphone  MM,  is 

j<d(t+AX/C  )  ju)(t-Ax/ C  ) 

PBK  =  A  e  0  +  A  e 

=  2A  cos(kAx) 

Therefore, 

P 

p^-  =  cos ( 2irf Ax/c  )  .  (4.8) 

kMM  0 

The  correction  factor,  in  decibels,  is  written  as 

B  =  20  log[cos(27rfAx/co]  .  (4.9) 

The  corrections  at  the  test  frequencies  are  as  follows:  -0.018  dB  at 
0.5  kHz;  -0.073  dB  at  1  kHz;  -0.293  dB  at  2  kHz;  and  -1.213  dB  at  4  kHz. 
The  adjusted  data  is  shown  in  Fig.  4.11.  The  average  sensitivity  is 
-108  dB  re  1  V/ybar,  a  value  that  is  within  1  dB  of  that  measured  by  the 
first  method.  The  validity  of  the  first  method,  therefore,  is  confirmed. 
We  used  the  sensitivity  determined  by  the  first  method  for  our  rise  time 
measurements. 

C.  Alignment  and  Procedure 

Alignment  of  the  spark  source  and  microphone  face  was  a 
critical  aspect  of  this  experiment.  The  better  the  alignment,  the  less 


I 


69 

critical  the  fact  that  our  microphone  is  not  a  true  point  receiver  (see 
the  discussion  of  microphone  directivity  in  the  section  on  the  receiving 
system).  The  alignment  was  done  with  the  light  beam  from  a  Hughes 
Model  31931H  laser.  First,  the  laser  beam  was  aligned  parallel  with  the 
axis  of  the  optical  bench.  This  was  done  by  shining  the  laser  beam  on  a 
screen  mounted  on  the  bench  and  adjusting  the  horizontal  position  of  the 
laser  until  the  beam  did  not  shift  when  the  screen  was  moved  down  the 
bench.  The  laser  was  then  secured  in  that  position.  The  spark  gap  and 
microphone  face  were  aligned  along  an  axis  formed  by  the  laser  beam.  The 
electrode  holder  and  laser  were  adjusted  vertically  until  the  beam  passed 
through  the  center  of  the  spark  gap  and  was  reflected  from  the  center  of 
the  microphone  face.  The  electrode  holder  was  then  adjusted,  by  eye, 
until  its  plane  was  perpendicular  to  the  laser  beam.  Finally  the  tilt  of 
the  microphone  face  was  adjusted  until  the  reflected  laser  beam  was 
superimposed  on  the  incident  beam.  We  were  then  sure  that  the  line 
between  the  center  of  the  spark  gap  and  microphone  face  was  perpendicular 
to  the  microphone  face. 

After  the  apparatus  was  aligned,  we  were  ready  to  take  data. 
First,  the  microphone-preamplifier  system  was  calibrated  (as  described 
above).  Next,  the  temperature  and  relative  humidity  were  measured  and 
recorded.  This  was  very  important  since  the  predictions  of  our  models 
depended  on  these  factors.  The  temperature,  as  measured  with  a  mercury 
thermometer,  remained  fairly  constant  all  year  long  at  75°  +2°F.  The 
humidity  was  measured  with  a  Brooklyn  direct  reading  hygrometer.  The 
range  of  relative  humidities  was  40-66%.  The  lower  humidities  generally 
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occurred  in  the  winter  months  and  the  higher  humidities  in  the  spring 
and  summer  months. 

Several  data  runs  were  then  taken.  Waveforms  were  captured  at 
several  source-receiver  distances.  The  same  averaging  procedure  used  in 
the  calibration  was  used  to  determine  which  waveforms  to  plot  on  the  X-Y 
recorder.  Some  waveforms  were  also  saved  on  a  disk.  The  waveforms  were 
then  analyzed  to  determine  the  rise  time.  Two  different  definitions  for 
rise  time  were  used  in  this  study,  as  shown  in  Fig.  4.12.  In  the  first 
definition,  the  rise  time  is  described  as  the  time  required  for  the  local 
pressure  to  increase  from  5%  to  95%  of  its  peak  value.  The  rise  time 
can  also  be  defined  in  terms  of  the  slope  of  the  initial  shock.  In  this 
case  the  rise  time  equals  the  amplitude  divided  by  the  slope  of  the  rise 
at  its  midpoint.  We  use  the  second  rise  time  definition  to  analyze  our 
data.  This  definition  is  chosen  because  it  is  the  one  most  currently  in 
use  in  the  literature. 

Some  sample  waveforms  are  shown  in  Fig.  4.13.  A  typical 
waveform  with  a  very  short  rise  time  is  shown  in  Fig.  4.13(a),  one  with 
a  much  longer  rise  time  in  Fig.  4.13(b). 

D.  Limits  on  Rise  Time  Measurement  Capability 

There  are  four  factors  which  limit  the  rise  time  measurement 
capability  of  our  measurement  system.  The  first  limitation  is  set  by  the 
oscilloscope.  Since  the  time  resolution  is  at  best  0.5  ysec/point,  we 
could  not  possibly  measure  rise  times  shorter  than  0.5  ysec.  The  second 
limitation  is  set  by  the  microphone.  The  shortest  rise  time  we  measured 
with  microphone  MM  is  about  0.65  ysec.  The  third  limitation  is  associated 
with  the  geometry  of  the  spark  source.  For  example,  the  electrode  holder 


(a)  RISE  TIME  DEFINITION  No.  1 


(b)  RISE  TIME  DEFINITION  No.  2 


FIGURE  4.12 

DEFINITION  OF  RISE  TIME 
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is  adjusted  perpendicular  to  the  laser  beam  by  eye.  If  the  gap  were 
slightly  tilted,  one  end  of  the  source  would  be  closer  to  the  microphone 
face  than  the  other.  If  the  trigger  electrode  is  not  centered  between 
the  electrodes  or  not  on  the  gap  axis,  there  may  be  a  path  length 
difference  from  various  parts  of  the  spark  to  the  microphone.  Lastly, 
the  electrodes  may  be  slightly  offset  from  each  other  Although  no  real 
analysis  of  the  effect  of  these  misalignments  was  made,  it  is  worthwhile 
noting  that  a  path  length  difference  corresponding  to  a  time  delay  (At) 
of  1.0  ysec  is 

AX  =  C  At 
0 

=  0.345  mm 

If  any  of  the  problems  discussed  above  results  in  a  path  length  difference 
of  0.345  mm,  our  rise  time  limitation  is  1.0  usee. 

Although  some  of  our  initial  waveforms  might  actually  have  rise 
times  smaller  than  1.0  usee,  most  of  our  data  are  not  affected  by  the 
limitations  imposed  by  the  equipment.  As  will  be  seen  in  Chapter  5,  most 
of  our  measured  rise  times  were  greater  than  1.0  ysec. 
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CHAPTER  5 

MEASUREMENTS  AND  COMPARISONS  WITH  THEORETICAL  PREDICTIONS 

A.  Introduction 

In  this  chapter  we  report  our  experimental  results  and  compare 
them  with  the  predictions  obtained  from  the  models  described  in  Chaps.  2 
and  3.  In  Section  B  our  results  are  compared  with  the  thermoviscous  and 
relaxation  model  predictions.  As  noted  in  Chap.  2  these  models  are  derived 
for  a  step  shock.  In  Section  C  we  present  a  comparison  of  the  waveforms 
computed  by  our  propagation  algorithm  with  those  from  our  experiments. 

Three  different  kinds  of  data  runs  were  taken.  First,  with  the 
spark  energy  and  gap  length  held  constant,  waveforms  were  captured  at 
specific  distances  between  10  and  550  cm.  In  the  second  kind  of  data  run 
we  varied  the  spark  energy,  gap  length,  and  propagation  distance  in  such 
a  way  that  the  half  duration  was  held  constant.  We  took  data  with  con¬ 
stant  ha^  duration,  because  we  wanted  to  check  whether  there  was  no 
dependence  of  rise  time  on  the  half  duration,  as  the  analytical  predic¬ 
tions  indicate.  In  the  last  kind  of  data  run  spark  energy,  gap  length, 
and  propagation  distance  were  varied  so  as  to  keep  the  peak  pressure 
constant. 

B.  Comparison  with  Thermoviscous  and  Relaxation  Models 

We  used  two  different  graphical  presentations  to  compare  the  data 
with  the  predictions  from  the  models.  The  first  is  a  plot  of  the  rise  time 
(tr)  versus  the  inverse  of  the  peak  pressure  (given  by  Eqs.  (2.7b)  and 
(2.14)).  This  method  was  chosen  because  the  rise  time  is  described  as  a 
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constant  times  the  inverse  of  the  peak  pressure  in  the  analytical  models. 

The  second  plot  is  of  the  slope  of  the  rise  at  its  midpoint  (m^)  versus  the 
peak  pressure  squared  (given  by  Eqs.  (2.7a)  and  (2.13a)).  We  chose  to  use 
this  method  after  we  compared  some  of  our  computed  waveforms  with  the  exper¬ 
imental  waveforms.  We  found  that  although  the  peak  pressures  of  the  two 
waveforms  were  not  exactly  the  same  (error  =8.5%) ,  the  slope  of  the  rises 
were.  Since  we  also  use  the  slope  of  the  rise  in  one  of  our  rise  time 
definitions,  this  method  seemed  reasonable. 

We  took  many  runs  with  both  constant  energy  and  constant  half 
duration.  The  runs  presented  here  are  representative  samples.  We  list 
the  different  runs  along  with  pertinent  information  about  them  in 
Tables  5.1(a)  and  5.1(b). 

1 .  Rise  Time  versus  Inverse  Peak  Pressure 
a.  Constant  Spark  Energy 

Basically,  we  found  that  the  rise  times  we  observed  are  closer  to 
those  predicted  by  the  thermoviscous  model  than  those  predicted  by  the 
relaxation  model.  Furthermore,  the  experimental  rise  times  agree  better 
with  TV2  than  TV1 .  A  plot  of  the  results  from  several  data  runs  is  shown 
in  Fig.  5.1.  As  noted  in  Chap.  4,  equipment  limitations  prevented  us  from 
measuring  a  rise  time  of  less  than  about  1.0  ysec.  This  accounts  for  the 
bunching  of  points  around  1.0  usee.  As  noted  in  Chap.  2,  the  predictions 
given  by  the  relaxation  model  are  only  valid  for  fully  dispersed  shocks, 
i.e.,  P<1  mbar.  Therefore,  the  relaxation  prediction  is  inapplicable  over 
most  of  the  plot  region.  However,  we  include  the  rise  time  prediction  for 
a  relaxing  medium  in  Figs.  5.1  and  5.2,  for  the  purpose  of  comparison  with 
the  thermoviscous  prediction. 
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TABLE  5.1 

SUMMARY  OF  DATA  RUNS 


(a)  CONSTANT  ENERGY 


SPARK  ENERGY 

GAP  1  FNGTH  TEMPERATURE 

RELATIVE  INITIAL 
HUMIDITY  PRESSURE 

INITIAL 

HALF 

DURATION 

MEASUREMENT 

RANGE 

(0) 

(mm) 

(°F) 

(*) 

(mbar) 

(nsec) 

(cm) 

0.061 

1 

73 

60 

1.03 

9.0 

50-300 

1.08 

- 

- 

- 

8.66 

7.5 

10-300 

1.96 

- 

75 

66 

7.33 

15.5 

30-325 

2.88 

- 

75 

66 

6.12 

15.5 

30-300 

2.88 

1 

75 

63-64 

5.39 

21.6 

30-550 

3.13 

2 

76 

60 

8.25 

17.0 

30-450 

4.81 

2 

73 

40 

4.50 

25.0 

60-450 

5.78 

2 

75 

60 

12.29 

24.5 

30-400 

6.48 

4 

77.5 

56 

14.41 

22.1 

30-450 

9.68 

1 

74 

63 

10.98 

19.7 

30-450 

(b)  CONSTANT 

HALF  DURATION 

HALF  DURATION 

MEASUREMENT 

RANGE 

RANGE 

SPARK  ENERGY 

GAP  LENGTH 
RANGE 

RELATIVE 
TEMPERATURE  HUMIDITY 

(usee) 

(cm) 

(J) 

(mm) 

(°F) 

(*> 

21.5 

50-400 

0.28-4.41 

- 

75 

66 

29.5 

75-550 

3.24-6. 

76 

- 

75 

66 

30.0 

50-450 

2,0-9.03 

2-3 

76 

52 

32.0 

100-500 

2.21-8. 

00 

- 

75 

66 

34.0 

100-450 

4.96-11 

.50 

2-4 

72 

50 

41.5 

75-600 

8.64-26 

.50 

3-5 

75 

66 

RISE  TIME  (t) 
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FIGURE  6.1 

FIRST  METHOD  OF  ANALYSIS: 
CONSTANT  SPARK  ENERGY 


b.  Constant  Half  Duration 
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Again  we  found  that  the  observed  rise  times  were  closer  to  those 
predicted  by  the  thermo viscous  model  (TV2 )  than  those  predicted  by  the 
relaxation  model.  We  captured  waveforms  with  half  durations  from  21.5 
to  41.5  /isec.  A  plot  of  the  results  from  three  constant  half  duration 
runs  is  shown  in  Fig.  5.2.  An  interesting  effect  that  may  be  seen  in  this 
plot  is  that  the  slope  of  the  curve  defined  by  the  data  points  seems  to 
increase  as  the  half  duration  increases.  This  is  rather  curious  because 
as  the  half  duration  increases  the  N  wave  begins  to  look  more  like  a 
step  shock.  One  conclusion  is  that  the  thermoviscous  model  is  not 
necessarily  better  than  the  relaxation  model.  It  may  be  that  the  half 
duration  is  not  long  enough  to  tell  whether  relaxation  has  an  effect  on 
the  shock  profile. 

c.  Constant  Amplitude 

Only  one  constant  amplitude  run  was  made.  We  captured  several 
0.62  mbar  peak  amplitude  waveforms  at  source-receiver  distances  of 
75-450  cm.  According  to  both  the  thermoviscous  and  relaxation  models, 
the  rise  times  for  all  of  these  waveforms  should  be  the  same.  We  did 
not  find  this  to  be  true.  We  measured  the  slope  m^  of  the  linear  decay 
portion  of  the  waveforms  and  plotted  the  rise  time  versus  slope;  see 
Fig.  5.3.  There  is  a  fairly  linear  relationship  between  the  rise  time 
and  the  linear  decay  slope;  as  the  slope  increases  the  rise  time 
decreases.  These  results  confirm  the  results  obtained  in  the  constant 
half  duration  runs.  As  the  half  duration  increases,  the  slope  of  the 
linear  decay  portion  decreases  and  the  rise  time  increases.  It  is 
possible  that  an  N  wave  as  long  as  a  sonic  boom  would  have  a  rise  time 
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equal  to  that  predicted  by  the  relaxation  model,  provided,  of  course, 
that  the  amplitude  were  low  enough. 

2.  Slope  of  the  Rise  (m-j )  versus  Peak  Pressure  Squared 

a.  Constant  Spark  Energy 

Again  we  find  that  our  data  fits  better  with  the  thermoviscous 
model  (TV2)  than  with  the  relaxation  model  (see  Fig.  5.4(a)).  In  this 
method  of  analysis  the  data  points  seem  to  fall  closer  together  and  the 
slope  of  a  least  squares  fit  is  similar  to  the  slope  of  the  thermoviscous 
lines.  There  are  a  few  points  at  very  high  peak  pressures  that  fall 
near  the  relaxation  line.  However,  the  relaxation  model  is  only  valid 
for  peak  pressures  under  1  mbar  (see  Chapter  2).  Therefore  we  cannot 
compare  these  points  with  the  results  from  the  relaxation  model.  A 
better  explanation  for  these  points  is  that  the  rise  is  so  fast  at  high 
amplitudes  that  we  are  approaching  the  limit  of  our  experimental  apparatus. 

b.  Constant  Half  Duration 

The  effects  of  a  lengthened  half  duration  are  not  as 

pronounced  in  this  method  of  analysis  as  in  the  previous  method.  A  plot 

of  the  same  constant  half  duration  runs  is  shown  in  Fig.  5.4(b).  However, 

it  is  still  evident  that  as  the  half  duration  increases,  the  data  seem  to 

fall  closer  to  the  relaxation  theory  line.  Again  we  have  those  few  data 

2 

points  above  3.2  mbar  that  seem  to  "bend"  towards  the  relaxation  line. 

This  is  probably  due  to  limitations  of  our  apparatus. 

C.  Comparison  with  Computer  Propagation  Algorithm 

The  results  from  th*  propagation  algorithm  can  only  be 
conveniently  compared  with  data  from  our  constant  spark  energy  runs.  We 
obtained  very  good  agreement  between  our  experimental  waveforms  and  the 
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computed  waveforms,  in  terms  not  only  of  general  wave  shape  but  of 
amplitude,  half  duration,  and  rise  time.  A  comparison  of  one  particular 
run  is  shown  in  Fig.  5.5.  If  we  plotted  the  data  from  the  computed  wave¬ 
forms  in  Fig.  5.1,  the  points  would  fall  very  close  to  the  x's.  Agreement 
is  poor  for  the  second  half  of  the  N  wave,  because  the  second  half  of  our 
experimental  N  waves  tend  to  be  contaminated  by  reflected  and/or  dif¬ 
fracted  signals.  The  maximum  error  in  amplitude  was  8.5%  and  the  maximum 
error  in  half  duration  was  2.0%.  Since  we  took  our  experimental  data 
using  several  sparks,  it  is  not  surprising  that  we  found  some  small  error. 
In  addition  to  spark  variability,  there  is  some  error  in  the  calibration 
process.  The  maximum  error  in  rise  time  is  27%.  This  is  a  relatively 
large  error.  Most  of  the  error  can  be  directly  attributed  to  error  in 
measuring  the  amplitude  because  the  slope  of  the  head  shock  m-j  was  found 
to  be  very  nearly  the  same  for  both  experimental  and  computed  waveforms. 
The  percentage  error  between  measured  and  predicted  rise  times  is  given  in 
Table  5.2  for  all  four  models. 


TABLE  5.2 

PERCENTAGE  DIFFERENCE  BETWEEN 
EXPERIMENTAL  AND  COMPUTED  RISE  TIMES 


r(cm) 

TV1 

TV2 

RELAXATION 

PROPAGATION 

ALGORITHM 

30 

81.91 

74.27 

-  155.6 

— 

52.84 

32.84 

-  566.45 

25.4 

200 

27.56 

-  3.10 

-  922.73 

26.6 

300 

18.61 

-15.87 

-1049.50 

23.4 

400 

15.66 

-23.92 

-1091 .43 

23.1 
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We  attribute  the  success  of  our  algorithm  to  the  fact  that  it  was 
able  to  accommodate  the  real  N  wave  signal.  A  step  shock  did  not  have  to 
be  assumed.  It  is  concluded  that  the  half  durations  of  the  N  waves  used 
in  our  experiment  are  not  long  enough  for  our  measurements  to  fairly  test 
the  step  shock  predictions. 

The  agreement  between  our  experimental  and  computed  waveforms  led 
us  to  suspect  that,  despite  the  apparently  poor  showing  in  Figs.  5.1  and 
5.2,  relaxation  does  have  an  appreciable  effect  on  the  waveforms.  In 
order  to  check  this  possibility,  we  ran  the  computer  algorithm  using  the 
following  absorptions:  thermoviscous  alone,  thermoviscous  and  relaxation 
(including  dispersion),  and  relaxation  (including  dispersion)  alone.  The 
computer  runs  were  made  using  the  same  input  waveform  and  initial  condi¬ 
tions  as  the  run  shown  in  Fig.  5.5.  Several  representative  waveforms  are 
shown  in  Fig.  5.6.  For  r=1.5  m,  the  slope  of  the  head  shock  for  the 
thermoviscous-only  and  relaxation- only  cases  is  steeper  than  the  slope  for 
the  combination  case.  However,  the  slope  for  the  thermoviscous-only  case 
is  much  closer  to  that  of  the  combination  case.  Therefore,  we  concluded 
that  thermoviscous  absorption  is  playing  a  dominant  role  at  r=1.5  m.  As 
the  propagation  distance  (r)  increases,  the  slope  of  the  head  shock  for 
the  relaxation-only  case  decreases  faster  than  that  for  the  thermoviscous- 
only  case.  By  the  time  we  reach  r=5.0  m  the  slope  of  the  head  shock  for 
the  relaxation-only  case  is  very  close  to  that  for  the  combination  case. 

In  the  thermoviscous-only  waveform,  however,  the  head  shock  is  much 
steeper.  We  therefore  conclude  that  relaxation  plays  the  dominant  role 
in  determining  the  head  shock  at.  r=5.0  m. 
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The  emerging  role  played  by  relaxation  as  the  N  wave  becomes 
weaker  can  also  be  seen  clearly  in  Fig.  5.7.  Here  we  plot  rise  time  ver¬ 
sus  propagation  distance  for  all  three  cases  and  for  the  experimental 
data.  As  the  propagation  distance  increases,  the  thermoviscous-only  rise 
time  increases  more  slowly,  while  the  relaxation-only  rise  time  increases 
more  quickly. 

An  explanation  for  our  findings  is  as  follows.  The  frequency 

spectrum  for  an  ideal  N  wave  (tr=0)  is  proportional  to  |jj(wT)|,  where  jj 

3Q 

is  the  spherical  Bessel  function  of  the  first  kind  and  order.  The  high 

frequency  envelope  of  the  spectrum  falls  off  at  6  dB/octave.  If  the  N 

wave  has  a  finite  rise  time,  however,  the  envelope  first  falls  at  6  dB/ 

octave  and  then  changes  to  a  12  dB/octave  slope.  The  slope  change  occurs 

40 

at  a  frequency  that  is  approximately  l/tr>  Thus,  as  the  rise  time  of 
the  N  wave  increases,  there  is  less  and  less  high  frequency  energy.  It 
is  well  known  that  for  frequencies  clearly  above  the  relaxation  frequen  y 
of  the  medium  the  dominant  absorption  mechanism  is  thermoviscous  whereas 
for  frequencies  clearly  below  the  relaxation  frequency  the  dominant  mech¬ 
anism  Is  relaxation.  Near  the  relaxation  frequency  there  is  a  sort  of 
"saddle"  region  in  which  both  mecha>  '  ire  important.  For  the  particu¬ 
lar  temperature  and  humidity  we  used  computer  run,  the  relaxation 

frequency  Is  about  74  kHz.  The  point  at  which  relaxation  absorption  and 
thermoviscous  absorption  are  equal  is  about  160  kHz,  which  corresponds  to 
a  rise  time  of  6.25  ysec.  Most  of  our  rise  times  fall  between  1.5  and 
6.25  usee.  At  each  end  one  of  the  absorption  mechanisms  dominate.  In 
the  middle,  both  mechanisms  seem  to  affect  the  rise  time  of  our  N  waves. 
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CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS 
FOR  FURTHER  WORK 

An  experiment  has  been  carried  out  to  measure  the  rise  time  of 
spherically  spreading  N  waves  produced  in  air  by  sparks.  The  main  ques¬ 
tion  is  whether  the  rise  time  was  determined  by  oxygen  relaxation.  For 
the  N  waves  measured— amplitudes  in  the  range  0.15-15  mbar  and  half 
durations  in  the  range  7.5-41.5  psec—a  conclusive  answer  to  this  question 
was  not  found.  It  did  appear,  however,  that  in  many  cases  the  amplitudes 
were  so  large  that  the  shocks  were,  at  best,  only  partly  dispersed. 

In  the  experiment,  N  waves  were  produced  by  means  of  an 
electric  spark  which  was  mounted  on  an  optical  bench.  Waveforms  were 
captured  at  specific  source-receiver  distances  with  a  condenser  micro¬ 
phone  (of  very  wide  bandwidth)  and  a  digital  oscilloscope.  Data  was 
taken  under  three  different  conditions:  constant  spark  energy,  constant 
half  duration,  and  constant  amplitude.  The  waveforms  captured  during 
these  runs  were  analyzed  to  determine  their  rise  time. 

We  have  reviewed  two  analytical  models  used  to  predict  the  rise 
time  of  a  step  shock,  the  thermoviscous  model  and  the  relaxation  model. 
Previous  investigators  have  used  these  models  to  predict  the  rise  time  of 
sonic  boom  N  waves.  Since  the  sonic  boom  has  a  very  long  half  duration 
(=0.2  sec),  a  step  shock  is  a  good  approximation  for  it.  Our  N  waves, 
however,  have  such  a  short  half  duration  (=20  psec)  that  a  step  shock  is 
hardly  a  good  approximation  for  them.  When  we  initially  compared  our 
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experimental  rise  times  with  those  predicted  by  the  models,  we  found  that 
the  experimental  data  seemed  to  fit  the  thermoviscous  model  better  than 
the  relaxation  model.  This  finding  should  not  be  weighed  heavily,  how¬ 
ever,  because  of  the  questionable  applicability  of  the  two  analytical 
model s . 

In  our  comparisons  using  the  constant  half  duration  runs  we 
noticed  a  curious  pattern.  As  the  half  duration  increased  so  did  the  rise 
time.  This  pattern  was  confirmed  when  we  obtained  a  similar  result  using 
the  constant  amplitude  run  data.  The  discovery  of  this  pattern  put  the 
initial  findings  described  above  in  question.  There  was  still  a  chance 
that  the  effects  of  oxygen  and/or  nitrogen  relaxation  were  important  but 
were  masked  by  the  assumption  that  an  N  wave  could  be  approximated  by  a 
step  shock. 

In  order  to  obtain  theoretical  predictions  valid  for  our  N 
waves,  we  took  a  direct  approach  to  the  problem.  We  used  a  wave  propaga¬ 
tion  computer  algorithm  in  which  the  effects  of  atmospheric  absorption 
(and  dispersion)  and  nonlinear  propagation  distortion  are  calculated 
directly.  This  algorithm  is  ideal  for  dealing  with  our  N  waves  because 
it  is  designed  to  accommodate  arbitrary  waveforms.  Given  an  experimen¬ 
tally  measured  N  wave  at  a  specified  distance,  our  wave  propagation 
algorithm  predicts  the  shape  of  the  N  wave  at  subsequent  distances. 
Atmospheric  absorption  was  calculated  using  a  digital  implementation  of 
the  American  National  Standard  Institute  (ANSI)  document  SI. 26-1978. 

The  absorption  calculation  Includes  the  effects  of  viscosity,  heat  conduc¬ 
tion,  and  vibrational  relaxation  of  oxygen  and  nitrogen  molecules.  Dis¬ 
persion  due  to  oxygen  relaxation  was  also  Included.  Nonlinear  distortion 
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was  calculated  using  standard  theory  of  finite  amplitude  waves.  The 
computed  waveforms  were  compared  with  the  experimental  waveforms.  The 
amplitude,  half  duration,  and  rise  time  were  in  good  agreement  for  all 
the  waveforms. 

It  is  concluded  that  the  existing  thermoviscous  and  relaxation 
models,  which  are  for  step  shocks,  cannot  be  used  to  predict  the  rise  time 
of  short  half-duration  N  waves.  As  for  oxygen  relaxation,  the  evidence  is 
that  it  does  have  some  effect  on  the  rise  time  of  our  N  waves.  Computer 
calculations  in  which  the  thermoviscous  and  relaxation  effects  could  be 
included  or  omitted  show  that  the  measured  rise  times  cannot  be  attributed 
to  either  mechanism  alone.  Both  mechanisms  are  important  for  the  N  waves 
for  our  experiments. 

Our  recommendation  for  further  work  is  twofold.  First,  true 
solutions  for  spherically  spreading  N  waves  should  be  obtained  for  the 
thermoviscous  and  relaxing  media.  Our  experimental  results  could  then  be 
used  to  test  the  solutions.  Second,  our  experiment  should  be  carried  out 
with  N  waves  of  longer  half  duration  and  smaller  amplitude.  A  greater 
range  of  results  will  be  useful  for  checking  theoretical  predictions. 


APPENDIX  A 


2  a 

DERIVATION  OF  THE  EQUATION  OF  STATE  FOR  A  RELAXING  MEDIUM 


A  sound  wave  passing  through  a  relaxing  medium  may  cause  a 
change  in  the  thermodynamic  equilibrium  of  the  fluid.  We  introduce  an 
"internal  coordinate"  5  to  account  for  these  changes.  If  |>=|>(p)  in  a 
nonreacting  medium,  then  for  a  relaxing  medium  we  must  write 

^=f(p,£)  .  (A.l) 

The  Taylor  series  expansion  of  Eq.  (A.l)  is 


P  -  (£-P0)  =  (|f)  p'  +  l 


, 2  . 
p  + 


(ll\ 

Wo 


(A.2) 


'5  VP  /5  '  ''P 

where  65=U-€0).  p,=(p-p0),  and  (a^/3p)?  is  the  square  of  the  frozen 
o 

sound  speed  (cf).  The  relaxation  (rate)  equation  is  usually  written 

_T  £  -  r 

T  dt  5  S> 

where  t  is  the  relaxation  time  of  the  medium  and  5Q  is  the  equilibrium 
state  of  the  internal  coordinate.  But  it  is  necessary  to  take  into 
account  the  dependence  of  50  on  density. 
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where  5  is  the  value  of  the  internal  coordinate  undisturbed  by  the 
00 

wave.  The  relaxation  equation  with  density  varying  because  of  the 
presence  of  the  sound  wave  is  therefore 


T  d£  _  .  ^0  , 

T  dt  5  5oo  dp  p  * 

or,  if  6£  stands  for  C-C00, 

«  '  -t5t  +  T?”'  •  <A-3) 

Substituting  Eq.  (A. 3}  into  the  last  term  on  the  right  hand  side  of 
Eq.  (A. 2),  we  find  that 


P  * 


(A. 4) 


If  €  is  constant,  i.e.,  the  fluid  is  in  equilibrium,  the  last  two  terms 
on  the  right  hand  side  of  Eq.  (A. 4)  become  zero.  Thus,  it  can  be  seen 
that 


is  equivalent  to  c. 


Substituting  this  result  Eq.  (A. 4)  now  becomes 


P  -  c^p 


(A. 5) 


Taking  the  time  derivative  of  Eq.  (A. 2)  and  multiplying  by  t  we  find  that 


T 


#  * 


1 

2 


(A. 6) 


Add  Eqs.  (A.5)  and  (A.6)  to  obtain 


1  +  T 


\  _  _  _2  i.l  /a2^\  .  d  \  ,2  ,  .2  do' 

;p-V  3t)p  TCf  dt  • 


Now,  add  and  subtract  (c£)  dp'/dt  from  the  right  hand  side  to  obtain 


+  T  h)p  s  {'  +  T  &)  cop' +  p" 


2  .  _  2  dp' 

I  +  mxc^  ji. 


o  ur  * 


2  2  2 

where  m*(Cf-c0)/c0.  This  is  Eq.  (2.7),  the  nonlinear  equation  of  state 
for  a  relaxing  but  inviscid,  non  heat-conducting  medium. 


,  »#  ' 


Sr, 


APPENDIX  B 


PROGRAM  AND  SUBROUTINE  LISTINGS 

Listings  for  Program  LB01 ,  Function  DABSORP  and  Subroutines 
RESAMP,  NWPLOT,  FQPLOT,  APWT,  and  TAYLOR  are  presented  in  this  appendix. 
Program  LBOI  is  used  to  compute  propagated  waveforms  from  a  given  direct 
waveform.  Function  DABSORP  is  used  to  compute  atmospheric  absorption 
and  dispersion  due  to  oxygen  relaxation  in  accordance  with  ANSI  Standard 
SI. 26-1978.  Subroutine  RESAMP  resamples  the  data  points  in  equal  time 
increments.  Subroutines  NWPLOT  and  FQPLOT  are  used  to  plot  the  computed 
waveforms  in  the  time  domain  and  the  frequency  domain,  respectively. 
Subroutines  APWT  and  TAYLOR  are  used  to  accurately  find  the  frequency 
components  of  the  N  wave. 
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